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EVALUATION 


The  work  presented  in  this  report  represents  an  effort  to  press  upon 
the  frontiers  of  knowledge  in  the  area  of  spectral  estimation,  spectral 
estimators,  and  the  performance  of  spectral  estimators.  To  date,  no 
researcher  has  developed  a  tractable  means  for  computing  the  performance  of 
a  maximum  entropy  spectral  estimator  in  the  presence  of  noise,  even  for  the 
type  of  noise  which  is  usually  most  convenient  mathematically:  Gaussian 
noise.  We  have  joined  the  ranks  of  these  researchers.  However,  we  have 
pressed  beyond  this,  point  of  disappointment  by  simulating  the  estimator 
against  noise,  signals,  and  interference,  and  by  composing  statistical 
analysis  of  the  estimator  outputs.  These  results  have  been  sufficiently 
encouraging  that  an  attempt  will  be  made  to  implement  the  maximum  entropy 
spectral  estimator  in  a  real  time  system,  and  to  validate  the  performance 
predicted  by  the  simulation. 

.  Vjo  Ji — 

KENNETH  E.  WILSON,  Capt,  USAF 
Project  Engineer 
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MAXIMUM  ENTROPY 


SPECTRAL  DEMODULATOR  INVESTIGATION 

1 .  INTRODUCTION  AND  OBJECTIVES 

A  transmitted  signal,  which  is  masked  by  noise  and,  possi¬ 
bly,  interference,  is  assumed  to  be  a  sinusoid  with  Hertz- 
frequency  which  varies  over  some  finite  set  F  of  positive  real 
numbers.  Let  "s(t)H  denote  the  value  of  the  resulting  continu¬ 
ous-time  signal  at  time  t.  Of  particular  interest  is  the  case 
where  #F  =  2;  that  is,  where  the  frequency  switches  back  and 
forth  between  two  distinct  values.  i 

The  object  of  this  study  is  to  evaluate  the  performance  of 
the  maximum  entropy  method  (=  MEM)  of  spectral  estimation  for 
short  segments  of  a  di screte-time  signal  which  results  from 
sampling  s(t)  at  uniformly  spaced  time  instants.  To  be  more 
specific,  it  is  desired  to  estimate  the  uncertainty  of  the  MEM 
estimates  of  the  transmitted  signal  frequencies,  by  obtaining 
confidence  intervals,  for  the  cases  (a)  received  signal  = 
transmitted  signal  +  noise  and  (b)  received  signal  =  transmit¬ 
ted  signal  +  noise  +  (purposeful)  interference. 

The  simple  case  (a)  with  F  =  {f^}  was  considered  almost 
exclusively.  The  (analytical)  problem  of  obtaining  confidence 
intervals  for  f^  ,  which  requires  determining  the  probability 
density  function  for  f-j  ,  appears  to  be  intractable.  Thus,  late 
in  the  program,  computer  simulation  was  used  to  study  the 
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sampling  variability  of  f  j .  The  results  of  this  simulation  are 
encouraging;  the  MEM  seems  to  perform  quite  well. 

2.  THE  MODEL 

In  general,  the  MEM  models  functional  values  as  the  output 
of  a  linear  discrete  system  (=  IDS);  that  is,  as  a  linear  coinbi 
nation  of  past  functional  values  (outputs  of  the  LDS)  and  past 
and  present  inputs  to  the  LOS.  This  leads  to  the  linear,  con¬ 
stant  coefficient,  difference  equation 

P  q 

(2.1)  f(kT)  «  -  l  a  f ( (k- j  )T)  +  G  l  b.u((k-i)T), 

j=l  J  i =0  1 

where  the  aj  and  b^  are  real  numbers,  b^  djF  l,  g  is  the  system 
gain  factor  (a  positive  real  number),  T  is  the  sampling  period 
(a  positive  real  number),  u((k-i)T)  Is  the  input  to  the  LDS  at 
time  (k-i)T,  and  k  is  any  positive  integer  such  that  the  func¬ 
tions  f  and  u  are  defined  at  the  indicated  times.  As  (2.1) 
enables  one  to  "predict"  f(mT)  from  f((m-l)T),  ...  ,  f((m-p)T), 
u(mT),  ...  ,  u((m-q)T),  the  name  "linear  prediction"  is  also 
associated  with  this  method. 

There  are  several  other  equivalent  representations  of  a 
LDS  in  addition  to  the  difference  equation  formulation  [1;  pp. 
85-86].  For  frequency-doma i n  considerations  ,  the  representa- 
ti  on 

S (z )  =  H  (z ) U ( z ) , 

where  S(z)  and  U(z)  are  the  z-transforms  of  s(kT)  and  u(kT) 
respectively  and  H(z)  is  a  rational  function  of  z  called  the 
"system  transfer  function,"  is  useful  [1;  pp.  220-282].  In 
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general,  H(z)  will  have  both  zeros  and  poles. 

We  model  a  received  signal  s  by  a  difference  equation  of 
the  form  (2.1)  with  q  =  0  and  G  =  1;  that  is,  we  assume 

P 

(2.2)  s(kT)  =-  J  a . s ( ( k- j )T )  +  u(kT) 

j=l  J 

for  all  appropriate  k.  For  the  model  (2.2),  it  can  be  shown 
that 


Thus  H  is  an  "all-pole"  transfer  function  with  p  poles,  namely, 

the  p  solutions  of  the  equation  1  +  £  p  a.z"J  =  0  or,  equiva- 

J  " !  J 

lently  (if  z  /  0,  as  is  required  by  the  definition  of  the  z- 
transform) , 

(2.3)  zP  +  a1zP"1  +  ...  +  ap_.|Z  +  ap  =  0. 

3.  ESTIMATION  OF  MODEL  PARAMETERS 

The  model  parameters  a^,  ...  ,  ap  in  (2.2)  are  approximated 
by  the  usual  type  of  least-squares  analysis  in  the  time-domain 
[2;  pp.  563-567].  For  a  given  positive  integer  m,  we  predict 
s(mT)  to  be 

_ _ „  P 

(3.1)  s(mT)  df-£  V  .s(  (m-j  )T) , 

j  =  l  J 

where  the  a",  are  chosen  so  as  to  minimize  an  appropriate  function 

J 

of  the  errors 

e ( kT )  df  s(kT)  -  C{M) 
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for  k  <  m.  In  the  case  of  a  deterministic  signal,  J  e2(kT)  is 
minimized  over  some  set  of  previous  samples.  In  case  we  mini¬ 
mize  the  sum  of  squared  errors 

m-1  N  - 1 

l  e2(kT)  (=  £  e2((m-N+k)T) 

k=m-N  k=0 


corresponding  to  the  preceding  N  samples,  the  technique  is  call¬ 
ed  the  "covariance  method"  of  linear  prediction  [2;  p.  564]. 

This  leads  to  the  system  of  linear  equations 


P 

(3.2)  l  a\ ( m , N ) $  ( m , N ) 

j=l  J  1J 


-*0i(m.N)  (i 


1,  2, 


P)  . 


where,  for  all  (i,j)  e  {1,  ...  ,  p}  x  { 1  ,  ...  ,  p } , 

N-1 

(3.3)  $..(nt,N)  df  l  s  ( (m-N+k-i  )T)s  ( (m-N  +  k-j  )T) , 

k  =  0 


for  determining  the  aj  which  yield  the  prediction  of  s(mT). 

From  (3.2)  and  (3.3),  we  see  that  the  N  +  p  consecutive  samples 
s((m-N-p)T),  ...  ,  s((m-l)T)  are  needed.  Thus  if  we  do  not 
allow  negative  arguments,  s((N+p+l)T)  is  the  first  sample  we 
can  predict. 

In  making  the  prediction  (3.1),  we  are  assuming  that  the 
input  u(mT)  is  completely  unknown,  which  is  often  the  case. 

4.  THE  ANALYTICAL  APPROACH 

The  analytical  determination  of  the  frequencies  f^  r  F 
exhibited  by  the  transmitted  signal  involves  (a)  calculating  the 
4>jj(m,N)  from  (3.3),  (b)  solving  the  system  of  linear  equations 
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(3.2)  for  the  IT,  ,  (c)  solving  the  polynomial  equation  (2.3) 

J 

with  "eTj"  in  place  of  "a^',  (d)  determining  e.  in  the  polar 
form  r^exp(±Je^)  of  each  of  the  complex  conjugate  pairs  of 
roots  of  (2.3),  and  (e)  multiplying  the  9^  by  an  appropriate 
real  number  to  obtain  the  f^ . 

As  indicated  in  Section  1,  an  attempt  was  made  to  handle 
analytically  the  simple  case  of  a  single  frequency  (F  =  { f ^ } ) 
signal  in  noise  with  no  interference;  that  is,  we  assume  the 
transmitted  signal  is 

(4.1)  Asin(2Tif^t) 
and  that  the  received  signal  is 

s(t)  =  A  sin(2nf.|t)  +  n(t), 

where  n(t)  is  a  zero-mean  Gaussian  noise  process  [3;  pp.  219- 
222]  which  is  uncorrelated  with  the  transmitted  signal  and  for 
which 

E[n(t)n(t+k)]  = 

Thus  the  samples  of  s  are  given  by 

(4.2)  s(kT)  =  A  si  n  (2irf  ^  kT )  +  n(kT). 

In  this  case,  a  2-pole  model  (p  =  2)  suffices. 

The  roots  of  (2.3)  with  p  =  2  and  "‘a."  in  place  of  "a." 

J  J 

are  z1  d_f  ( -a^  +  /~J)/2  and  z-j  d_f  ( -a^  -  /“d)/2,  where  d  df  2 
Aa^  It  can  be  shown  that  the  radian  argument  in  the  polar  form 
of  z-j  is  the  radian-frequency  f^  in  (4.2).  (In  fact,  this  poten 


fa2  (k  =  0), 

\  0  (k  f  0). 


1 1  a  1  for  ferreting  out  the  frequencies  of  a  sinusoidal  signal 
Is  the  reason  linear  prediction  Is  useful  to  us,  rather  than 
Its  predictive  powers.  We  can  always  "wait"  T  seconds  and 
measure  s(mT).)  Thus,  as  f  >  0,  z^  must  be  a  non- real ( compl ex ) 


number;  that  is. 

we  must  have  a^2  -  4a^  <  0. 

Hence  we  have 

z 

!  =  -  (a^/2)  +  J  (/  4r2  -  t}2)/2 

• 

where  4a^  -  a~^  2 

>  0.  It  is  easy  to  show  that 

z^  has  polar  form 

r^exp(Jt)^  ) ,  where 

(4.3) 

",  -  /r2 

and 

r  */2 

(^  =  0), 

(4.4)  e1  =  < 

tan  ^  ( -  /  4  &  2  ^  2  /  a^  )  +  n 

( a"i  >  0 ) , 

1 

tar>-,(-/  4a,  *  T.  s/a"  ) 

^  2)1 

O 

V 

1* 

If  we  measure  frequencies  in  Hertz  and  require 

then 

that  F  C  (0,6000) 

(4.5) 

f]  =  ( 6000/ it  ) 0 1  . 

Now,  'aj  and 

a^  in  the  expression  for  f  can  be  obtained 

explicitly  by  solving  the  system  of  equations 

(3.2)  with  p  =  2. 

By  Cramer's  Rule 

,  we  have 

(4.6) 

and 

al  =  ^*12*02  '  *22*01  )/A 

(4.7) 

a2  =  (*21*01  "  *02*11 )/A’ 

I 
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where 


(4.8)  A  df  ^ 1 1 ^ 2  ~  ^12^21 

and  the  are  given  by  (3.3)  and  involve  the  s(kT)  as  given  by 
(4.2).  Even  in  this  very  special  case,  the  complicated  chain  of 
operations  linking  the  assumptions  about  the  transmitted  signal 
and  the  noise  with  the  frequency  f^  makes  it  difficult  to  deter¬ 
mine  the  probability  density  function  for  f^. 

In  case  *F  *  2,  a  4-pole  model  is  appropriate  and  the  equa¬ 
tion  (2.3)  is  quartic.  The  four  roots  of  a  quartic  equation  can 
be  given  explicitly  in  terms  of  radicals  and  the  coefficients 
a.j,  ...  ,  a^;  however,  the  expressions  for  these  roots  are 
extremely  complicated.  Of  course,  if  #F  >  2,  a  model  with  at 
least  six  poles  is  required  and  no  analog  of  the  quadratic  and 
quartic  formulas  exists  for  the  equation  (2.3)  in  such  cases. 
Thus,  for  these  cases,  we  can  not  mimic  the  treatment  of  the 
2-pole  case. 


5.  TH£  COMPUTER  SIMULATION 


The  first  part  of  the  program  involved  an  extensive  inves¬ 
tigation  of  linear  prediction  techniques  and  their  application 
to  communications  theory  [1].  This  investigation  Included  the 
analytical  effort  described  in  Section  4.  The  second  part  of 
the  program  Involved  the  development  of  a  computer  simulation 
to  determine  the  sampling  variability  of  the  f^. 
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In  Sections  1  and  4,  we  did  not  elaborate  on  the  term 
"noise."  Noise  is  a  fundamental  limitation  on  the  performance 
of  physical  systems  such  as  radar  devices.  Some  possible 
sources  of  noise  are  clutter,  cosmic  radiation,  and  thermal 
motion  of  the  electrons  and  ions  in  the  receiver  components 
and  the  antenna  surroundings  [3;  pp.  3-5]. 

There  are  other  limitations  on  performance  in  signal¬ 
processing  that  are  based  on  the  fact  that  the  values  of  vari¬ 
ables  in  models  of  real-world  systems  are,  typically,  real 
numbers  with  decimal  representations  which  are  non-termina¬ 
ting  or  terminate  only  after  a  large  number  of  digits,  whereas, 
digital  si gnal -processi ng  equipment  (which  for  various  reasons 
has  largely  supplanted  analog  equipment)  seldom  allows  repre¬ 
sentations  of  numbers  (or  other  symbols)  with  more  than  64 
precision  bits.  One  such  limitation,  "quantization  noise," 
results  from  use  of  an  anal og-to-di gi tal  converter  (=  ADC)  in 
sampling  the  continuous-time  signal.  Another,  "roundoff  (or 
chopping)  noise,"  results  from  rounding  (or  chopping)  sums  and 
products  to  fit  the  computer's  word  length. 

In  the  simulating  done  to  date,  an  effort  has  been  made  to 
assess  the  sampling  variability  due  to  quantization  noise  and 
roundoff  noise.  Such  noise  is  always  present  as  we  can  not  do 
infinite-precision  sampling  and  arithmetic.  Further  work  is 
needed  to  determine  the  effect  of  purposeful  man-made  interfer¬ 
ence  and  the  type  of  noise  cited  in  the  second  paragraph  of  this 
section.  Also,  we  have  considered  only  the  2-pole  model  dis¬ 
cussed  in  Section  4. 
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Some  features  of  the  simulation  which  look  peculalr  owe  to 
the  fact  that,  previously,  some  processing  of  actual  signals  was 
done  at  RADC.  An  8 - b 1 1  ADC  was  used  to  produce  8 - b  1 1  approxima¬ 
tions  of  signal  values  and  blocks  of  2048  of  these  were  stored 
on  magnetic  tape.  Later,  this  Information  was  computer-process¬ 
ed.  The  programs  used  in  the  present  study  include  portions  of 
^ the  previously-used  program. 

Below  is  a  list  of  key  variables  in  the  simulation  programs 
along  with  the  corresponding  variable  (if  any)  in  Section  4  and 
its  interpretation: 


TSA  (T) 

sampling  period  (1/12000  sec.); 

FI  (f, ) 

Hertz-frequency  of  the  transmitted  signal; 

A  (A) 

amplitude  of  the  transmitted  signal 

NP  (p) 

number  of  poles; 

N  (N) 

number  of  error  terms  in  the  minimization 
process  for  determining  the  prediction 
coefficients  in  (3.1); 

NN  (=  N  +  p) 

the  number  of  previous  samples  needed  to 
predict  s(mT) 

NB  (  ) 

number  of  ADC  bits; 

S(J)  ( s ( jT  ) ) 

value  of  the  received  signal  at  time  jT ; 

P(I,0)  (♦1  (m,N)) 

coefficient  of  a^  in  equation  i  of  the 
linear  system  *  Is  defined  In  (3.3); 

P0(I)  U0i(m,N)) 

constant  on  the  right  in  equation  1  of 
the  linear  system  (3.2),  $  Is  defined  In  (3 

R(K)  ( r] ) 

magnitude  of  the  root  of  (2.3)  in  the 
upper  half-plane; 

F(K)  (f, ) 

Hertz- frequency  corresponding  to  the  root 
of  (2.3)  In  the  upper  half-plane; 

FC  (  ) 

the  radian-Hertz  conversion  factor  6000/n 
in  (4.5); 

DEL  (a) 

the  determinant  of  the  system  (3.3)  (see 
(4.8)); 

A1  [T}) 

coefficient  in  equation  (2.3)  with  p  =  2; 

A2  ( a 2 )  coefficient  in  equation  (2.3)  with  p  =  2; 

AM  arithmetic  mean  of  F  df  { f ( mT ) :  m  e 

{6,  7 .  2048}}; 

SO  standard  deviation  of  F,  the  frequencies 


We  assume  henceforth  that  the  amplitude  A  (measured  in  volts) 

of  the  transmitted  signal  is  in  [-5,  5].  In  the  simulation  of 

the  ADCs  behavior,  the  signal  S ( J )  in  the  interval  [-5,  5],  of 

length  10,  is  mapped  into  the  integer  Interval  [-2NB_1,  2NB_1] 

N  B 

by  following  xi — *  (2  /10)x  by  chopping  to  an  integer,  and  the 

integer  is  mapped  back  into  a  computer  real  number  in  [-5,  5]  by 
x  * — *(10/2NB)x. 

In  a  certain  average  sense,  binary  representations  of  numbers 
have  log^flO)  (*  3.32)  times  as  many  symbols  as  do  the  correspond¬ 
ing  decimal  representations  when  both  representations  terminate. 
Thus  it  was  apparent  to  the  writer  that  the  8-bit  ADC  was  not 
adequate  for  dealing  with  frequencies  in  (0,  6000). 

One  of  the  programs  used  in  the  simulation  appears  in  the 
Appendix.  Its  purpose  is  to  determine  the  effect  of  the  number 
of  ADC  bits  on  the  performance  of  the  model.  The  results,  given 
in  Table  1,  of  executing  this  program  confirm  the  conclusion 
about  the  inadequacy  of  the  8- b 1 t  ADC. 
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This  program  uses  5  samples  (N=3,  p=2  )  to  form  each  estimate 


of  the  frequency  (f,  ).  The  frequency  is  estimated  2043  times, 
using  2048  samples.  These  samples  are  generated  with  a  specified 
number  of  bits  of  ADC  precision  (NB)  using  the  procedure  outlined 
previously.  The  mean  (AM)  and  standard  deviation  (SO)  are 
computed  in  subroutine  STATS  from  the  2043  frequency  estimates 
using  the  following  formulas: 

,  2043 

AM  =  2043  f i 

,  2043  2 

SD  =  twt  (fi  -  3250  } 

Table  2  shows  the  variation  of  AM  and  SD  with  N,  the  number 
of  error  terms  used  in  the  minimization.  In  the  case  NB  =  14, 
central  processor  (=  CPU)  time  is  given  also.  These  times  include 
the  execution  time  for  the  statistical  calculations.  Thus  absolute 
differences  are  meaningful  to  the  determination  of  the  effect  on 
execution  time  of  changing  the  value  of  N,  whereas  relative 


TABLE  1 


Mean  AM  and  standard  deviation  SD  of  predictions  for 
the  number  of  ADC  bits  NB  =  2,  3,  ...  ,  16,  27*  with 
A*  1.5,  N  3  3,  and  FI  3  3250. 


NB  AM  SD 


2 

6000.0000 

3 

3497.7974 

4 

3252.4888 

5 

3248.9677 

6 

3249.5983 

7 

3249.6557 

8 

3249.9485 

9 

3249.9887 

10 

3249.9991 

11 

3250.0001 

12 

3250.0004 

13 

3249.9985 

14 

3249.9960 

15 

3250.0002 

16 

3250.0008 

27 

3250.0001 

2749.9994 
1677.3801 
226.23133 
78.865247 
33.414235 
21 .808373 
1 1 .663445 
4.8569507 
2.4628091 
0.77203118 
0.22295413 
0.18255417 
0.15308648 
0.075419844 
0.034125918 
0.0066758151 


differences  are  not.  It  appears  that  adding  1  to  the  value  of  N 
increases  execution  time  by  approximately  (1/3) (0.0001  )  hr.  or 
0.12  sec.  Note  that  execution  of  the  program  Involves  2043 
(=  2048  -  (3  +  2))  calculations  of  an  f ^ ,  Thus  the  added  time 
for  calculation  of  each  value  of  f^  is  approximately  5-10“^  sec. 
Execution  times  are  of  interest  as  the  ultimate  goal  is  real¬ 
time  Implementation  of  the  model. 

♦Here  we  have  used  27  precision-bit  floating-point  representations 
and  arithmetic  with  the  ADC  simulation  portion  of  the  program 
deleted. 


TABLE  2 


Mean  AM  and  standard  deviation  SO 
the  number  of  error  terms  used  in 
=2,3,...  ,  12  with  A  =  1 .5,  FI 
10,  12,  14,  16,  and  27. 


of  predictions  for 
the  minimization  N 
=  3250,  and  NB  =  8, 


NB  =  8 


N 

AM 

SD 

2 

3250.0737 

15.632812 

3 

3249.9485 

1 1.663445 

4 

3249.9371 

9.3071191 

5 

3249.8821 

7.1128755 

6 

3249.8813 

6.0597079 

7 

3249.8710 

5.2067635 

8 

3249.8677 

4.4239790 

9 

3249.8625 

3.3021908 

10 

3249.8594 

2.8550067 

11 

3249.8582 

2.2166742 

12 

3249.8565 

2.0223015 

NB  =  10 


N 

AM 

SD 

2 

3250.0018 

2.9927573 

3 

3249.9991 

2.4628091 

4 

3249.9970 

2.0885355 

5 

3249.9975 

1.7026238 

6 

3249.9979 

1.6084557 

7 

3249.9978 

1.2938950 

8 

3249.9922 

1.0680774 

9 

3249.9925 

0.73630817 

10 

3249.9945 

0.56215978 

11 

3249.9921 

0.39859506 

12 

3249.9906 

0.35790581 

—-ovooo^io^cn^ojro 


NB  =  12 


3249.9978  0.28106676 

3250.0004  0.22295413 

3250.0015  0.14402156 

3250.0015  0.13972569 

3250.0001  0.081098709 

3250.0009  0.067171544 

3249.9994  0.059041822 

3250.0011  0.044649824 

3250.0009  0.047537731 

3249.9999  0.052555363 

3250.0009  0.046175324 


CPU  TIME 


3249.9988  0.21686259  0.0009 

3249.9960  0.15309648  0.0009 

3249.9997  0.11826708  0.0010 

3249.9983  0.10109097  0.0010 

3249.9958  0.066045111  0.0010 

3249.9998  0.059400090  0.0011 

3249.9980  0.052958412  0.0011 

3249.9972  0.032415771  0.0011 

3250.0021  0.037008012  0.0012 

3249.9973  0.038732876  0.0012 

3250.0007  0.033104122  0.0012 


j  *  - ■» 


NB  =  16 


N 

AM 

SO 

2 

3249.9999 

0.045065880 

3 

3250.0008 

0.034125918 

4 

3249.9995 

0.024391433 

5 

3249.9991 

0.021149200 

6 

3250.0010 

0.019371934 

7 

3249.9993 

0.015833400 

8 

3249.9985 

0.013296252 

9 

3249.9994 

0.011714947 

10 

3249.9999 

0.010651756 

11 

3249.9995 

0.0083492643 

12 

3250.0000 

0.0071696392 

NB  =  27 


N 

AM 

SO 

2 

3250.0003 

0.0087628514 

3 

3250.0001 

0.0066758151 

4 

3250.0000 

0.0042304078 

5 

3250.0000 

0.0027758344 

6 

3250.0000 

0.0017828080 

7 

3250.0000 

0.0020461476 

8 

3250.0000 

0.0021599298 

9 

3250.0000 

0.0019517387 

10 

3250.0000 

0.0015333511 

11 

3250.0000 

0.00091071260 

12 

3250.0000 

0.0010223139 
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The  next  table  shows  the  variation  of  AM,  SD,  and  the 
si gnal -to- (quanti zati on  and  roundoff)  noise  ratio  SNR,  given 
by  20* log10(2NBA/10) ,  with  A. 


TABLE  3 


Mean  AM,  standard  deviation  SD,  and  si  gnal -to- (quan¬ 
tization  and  roundoff)  noise  ratio  SNR  for  the  ampli¬ 
tude  of  the  transmitted  signal  A  =  0.5,  1.0,  1.5,  ... 
5.0  with  N  =  3,  FI  =  3250,  and  NB  =  10,  12,  14,  and  16. 


NB  =  10 


A 

AM 

SD 

SNR 

0.5 

3249.9668 

6.9402777 

34.185400 

1.0 

3249.9982 

2.0219451 

40.205999 

1.5 

3249.9991 

2.4628091 

43.727824 

2.0 

3249.9991 

1.5623139 

46.226599 

2.5 

3249.9962 

1 .5531899 

48.164799 

3.0 

3250.0001 

0.77203118 

49.748425 

3.5 

3250. 0002 

1.3622908 

51 .087360 

4.0 

3250.0012 

0.72204416 

52.247199 

4.5 

3250.0019 

0.42808743 

53.270249 

5.0 

3249.9961 

0.98476960 

54.185400 

NB  =  12 

A 

AM 

SD 

SNR 

0.5 

3249.9991 

1.5623139 

46.226599 

1.0 

3250.0012 

0.72204416 

52.247199 

1  .  5 

3250.0004 

0.22295413 

55.769024 

2.0 

3249.9994 

0.53745331 

58.267799 

2.5 

3249.9979 

0.44305741 

60.205999 

3.0 

3249.9985 

0.18255417 

61 .789624 

3.5 

3250.0007 

0.20943003 

63.128560 

4.0 

3249.9981 

0.27760254 

64.288399 

4.5 

3250.0016 

0.16504455 

65.311449 

5.0 

3250.0001 

0.16447618 

66.226500 
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NB  =  14 


A 

AM 

SD 

SNR 

0.5 

3249.9994 

0.53745331 

58. 267799 

1.0 

3249.9981 

0.27760254 

64.288399 

1  .  5 

3249.9960 

0.15309648 

67.810224 

2.0 

3250.0001 

0.12638791 

70.308999 

2.5 

3250.0036 

0.086964880 

72.247199 

3.0 

3250.0002 

0.075419844 

73.830824 

3.5 

3249.9987 

0.056363066 

75.169760 

4.0 

3250.0002 

0.056470926 

76.329598 

4.5 

3250.0023 

0.025813821 

77.352649 

5.0 

3249.9998 

0.049159881 

78.267799 

NB  =  16 


A 

AM 

SD 

SNR 

0.5 

3250.0001 

0.12638791 

70.308999 

1.0 

3250.0002 

0.056470926 

76.329598 

1  .  5 

3250.0008 

0.034125918 

79.851424 

2.0 

3249.9998 

0.033097830 

82.350199 

2.5 

3250.0016 

0.024226625 

84.288399 

3.0 

3250.0005 

0.017549278 

85.872024 

3.5 

3249.9992 

0.021185211 

87.210959 

4.0 

3249.9992 

0.018053080 

88.370798 

4.5 

3250.0004 

0.015683560 

89.393849 

5.0 

3249.9998 

0.015091892 

90.308999 

Table  4  shows  the  variation  of  AM  and  SO  with  the  frequency 
FI  of  the  transmitted  signal.  As  there  are  14-bit  ADCs  avail¬ 
able  commercially  and  model  performance  is  quite  good  with  NB  = 
14  according  to  Table  1,  we  let  NB  =  14. 

The  strange  variation  of  SD  with  FI  results  from  the  sam¬ 
pling  process  and  finite  precision  arithmetic.  As  the  possible 
signal  frequencies  are  less  than  6000  Hz.,  we  must  sample  at  the 
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TABLE  4 


Mean  AM  and  standard  deviation  SD  of  predictions  for 

FI  =  1  25  ,  250  ,  375  .  5875  with  NB  =  14,  N  =  3, 

and  A  =  1.5. 


FI 

AM 

SD 

125 

119.24084 

36.594588 

250 

249.51668 

13.343529 

375 

374.99050 

2.1931926 

500 

499.99242 

2.2354304 

625 

624.99677 

1.6018316 

750 

750.00167 

0.41421278 

875 

874.99947 

0.93162946 

1000 

999.99837 

0.73667024 

1125 

1124.9997 

0.35508527 

1250 

1250.0013 

0.36510814 

1375 

1374.9998 

0.33071683 

1500 

1499.9985 

0.24333926 

1625 

1624.9985 

0.24277011 

1  750 

1750.0004 

0. 1  09731  79 

1875 

1875.0006 

0.10363140 

2000 

2000.0000 

0. 

2125 

2124.9995 

0.18358279 

2250 

2250.0019 

0.073796298 

2375 

2375.0001 

0.21530845 

2500 

2499.9964 

0.090456136 

2625 

2625.0006 

0.12147366 

2750 

2750.0028 

0.15310039 

2875 

2875.0002 

0.1831 3020 

3000 

3000.0000 

0. 

3125 

3124.9978 

0.18266767 

3250 

3249.9960 

0.15309648 

3375 

3374.9995 

0.12147751 

3500 

3500.0056 

0.090457460 

3625 

3625.0008 

0.21415340 

3750 

3749.9958 

0.073787051 

3875 

3874.9984 

0.18436528 

4000 

4000.0000 

0. 

4125 

4124.9997 

0.10363143 

4250 

4249.9976 

0.10972426 

4375 

4374.9999 

0.24238029 

4500 

4499.9988 

0.24332672 

4625 

4624.9993 

0.33069 280 

4750 

4750.0037 

0.36511440 

4875 

4875.0005 

0.35509441 
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FI 

AM 

SD 

5000 

4999.9958 

0.73666593 

5125 

51 24.9982 

0.90562031 

5250 

5250.0022 

0.41419801 

5375 

5375.0008 

1  .5328385 

5500 

5500.0052 

2.2354247 

5625 

5625.0036 

2.1931942 

5750 

5750.4873 

13.343526 

5875 

5880.7836 

36.559757 

rate  of  at  least  12000  samples  per  second  according  to  the  Sam¬ 
pling  Theorem  [1;  p.  291].  Thus  we  are  led  to  letting  T  = 
1/12000.  For  FI  =  2000  (respectively,  3000,  4000),  we  have  a  df 
2kirFl/12000  =  k tt / 3  (respectively,  kir/2,  2kn/3)  and  sin(a)  =  .50 
(respectively,  1.0,  .50)  in  exact  arithmetic.  As  the  library 
SIN  function  is  accurate  to  8  decimal  digits,  sin(a)  is  exact 
and  1.5sin(a)  is  also  exact  with  two  significant  decimal  digits. 
In  these  cases,  it  can  be  shown  that  the  first  step  in  the  simu¬ 
lation  of  the  ADCs  behavior  leads  to  an  integer  if  NB  =  14;  thus 
the  ADC  simulation  leads  to  the  exact  frequency  FI  as  the  pre¬ 
dicted  frequency.  Thus,  the  value  0  for  SD  in  case  FI  =  2000, 
3000,  or  4000  is  explained.  Similar  reasons  for  other  irregu¬ 
larities  can  be  given. 

The  program  leading  to  Table  4  was  executed  with  FI  ranging 
from  125  to  5875  by  steps  of  25;  however,  for  obvious  reasons 
not  all  of  these  values  appear  in  the  table.  The  original  at¬ 
tempt  at  execution  for  FI  ranging  from  25  to  5975  lead  to  execu¬ 
tion-time  errors.  Further  work  showed  that  the  trouble  begins 
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somewhere  In  the  intervals  [100,  125)  and  (5875,  6000).  Obvi¬ 
ously,  as  FI  approaches  0f  (6000'),  the  radian  measure  of 
in  (4.4)  must  approach  0+  («*).  Thus  the  discriminant  d  df 
4a^  -  a^2  must  approach  0  in  both  cases.  Printing  values  of 
some  intermediate  variables  for  the  bad  cases  showed  that  d 
took  on  negative  values  with  increasing  frequency  as  FI  got 
closer  and  closer  to  0+  (6000').  Quantization  and  roundoff 
noise  had  taken  its  toll. 

The  next  table  shows  the  variation  of  AM  and  SD  with  block 
size.  Ultimately,  a  block  size  smaller  than  2048  may  be  used. 
We  note  that  the  performance  of  the  model  is  reasonably  uni¬ 
form  over  block  size. 


TABLE  5 


Mean  AM  and  standard  deviation  SD  of  predictions  for 
block  size  =  8,  16,  32,  64,  128,  256,  512,  1024,  and 
2048  with  NB  =  14,  FI  =  3250,  and  A  =  1.5. 


AM 


SD 


BLOCK  SIZE 


8 

16 

32 

64 

128 

256 

512 

1024 

2048 


3249.9163 

3249.9741 

3249.9907 

3249.9952 

3249.9980 

3249.9989 

3250.0002 

3249.9998 

3249.9960 


0. 10430552 
0.1 9046609 
0.14852447 
0.16077075 
0.15215203 
0.15498016 
0.15291365 
0.1 5360796 
0.15309648 


Finally,  relative  frequencies  and  cumulative  relative  fre¬ 
quencies  of  the  predictions  were  obtained  in  the  absence  of  the 
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probability  density  function  (=  PDF)  and  the  distribution 
function  (=  DF). 


TABLE  6 


Approximations  of  the  probability 
PDF  and  the  distribution  function 
tions  with  NB  =  14,  N  =  3,  A  =  1.5 

density  function 
OF  of  the  predic 
,  and  FI  =  3250. 

PREDICTION 

APPROX.  PDF 

APPROX.  DF 

INTERVAL 

ON  INTERVAL 

AT  RIGHT 

3249.45-3249.50 

0.041605482 

0.041605482 

3249.50-3249.55 

0. 

0.041605482 

3249.55-3249.60 

0. 

0.041605482 

3249.60-3249.65 

0. 

0.041605482 

3249.65-3249.70 

0. 

0.041605482 

3249.70-3249.75 

0. 

0.041605482 

3249.75-3249.80 

0.16691140 

0.20851689 

3249.80-3249.85 

0. 

0.20851689 

3249.85-3249.90 

0.083700441 

0.29221733 

3249.90-3249.95 

0.16691140 

0.45912873 

3249.95-3250.00 

0.20802741 

0.66715614 

3250.00-3250.05 

0. 

0.66715614 

3250.05-3250. 10 

0.16642193 

0.83357808 

3250.10-3250.15 

0.083210965 

0.91078904 

3250.15-3250.20 

0.083210965 

1 .00000000 

The  "grainy"  nature  of  the  above  results  is  again  the  effect 
of  quantization  and  roundoff.  To  partially  negate  these  effects 
the  program  was  run  again  with  the  quantization  portion  missing 
so  that  27-bit  f  1  oati ng-poi nt  representations  and  arithmetic 
prevail  . 
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TABLE  7 


Approximations  of  the  probability  density  function 
PDF  and  the  distribution  function  DF  of  the  predic¬ 
tions  with  N  =  3,  A  =  1.5,  FI  =  3250,  and  27-bit 
floating-point  representations  and  arithmetic. 


PREDICTION 

APPROX.  PDF 

APPROX.  DF 

INTERVAL 

ON  INTERVAL 

AT  RIGHT 

3249.985-3249.990 

0.12922173 

0.12922173 

3249.990-3249.995 

0.0P1742535 

0.21096427 

3249.995-3250.000 

0. 21*095937 

0.49192364 

3250.000-3250.005 

0.36^14929 

0.85707294 

3250.005-3250.010 

0.064121390 

0.92119432 

3250.010-3250.015 

0.042094959 

0.95790504 

3250.015-3250.020 

0.042094959 

1 .00000000 

6.  CONCLUSIONS  AND  RECOMMENDATIONS 

As  indicated  earlier,  the  problem  of  obtaining  confidence 
intervals  for  the  f^  e  F  appears  to  be  intractable,  even  for  the 
simple  case  of  a  single-frequency  signal  (F  =  {f-j})  in  noise. 
However,  simulation  of  the  s i ngl e- frequency  case,  with  the 
only  noise  being  quantization  and  roundoff  noise,  has  shown 
that  the  MEM  performs  quite  well  in  terms  of  sampling  variabil¬ 
ity  of  f  -j  . 

If  results  similar  to  those  in  Table  1  hold  for  the  two- 
frequency  signal  simulation  and  the  frequencies  differ  by,  say, 
50  Hz.,  it  is  obvious  (as  SD  =  11.7)  that  an  8-bit  ADC  is  inade¬ 
quate  and  that  an  ADC  giving  12-16  bits  is  desirable.  A  catalog 
search  has  shown  that  12-bit  ADCs  with  2usec.  sampling  time  and 
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and  1 4- b 1 t  ADCs  with  50  ysec.  sampling  time  are  available  com¬ 
mercially. 

From  Table  2  with  NB  *  14,  we  see  that  model  behavior  Is 
good  with  N  *  3  but  that  SD  can  be  decreased  by  a  factor  of 
approximately  4/5  (that  is,  to  1/5  of  its  value  for  N  *  3)  by 
increasing  N  to  12  at  the  cost  of  Increasing  CPU  time  by  a  fac¬ 
tor  of  approximately  1/3.  If  the  s 1 ngl e- frequency  case  were 
of  practical  Interest,  It  would  not  seem  desirable  to  Increase 
N  in  view  of  the  desire  to  operate  In  real-time  and  the  good 
performance  of  the  model  with  N  =  3  (and  NB  =  12-14). 

From  Table  3  with  NB  =  12  (or  14),  we  see  that  model  per¬ 
formance  is  good  with  transmitted  signal  amplitude  A  =  1.5. 
However,  SD  can  be  decreased  by  a  factor  of  approximately  3/4 
(or  5/6)  by  increasing  A  to  4.5. 

From  Table  4,  we  see  that  model  performance  is  quite  good 
for  fj  e  (1000,  5000)  but  deteriorates  rapidly  as  fj  approaches 
0+  and  6000".  Again,  if  similar  results  hold  for  the  two-fre¬ 
quency  signal  case,  values  of  f^  and  f^  in  [1500,  4500]  would 
give  excellent  results. 

Table  5  shows  that  model  performance  is  quite  uniform  and 
good  for  block  sizes  ranging  from  128  to  2048.  Thus  a  reduction 
of  block  size  to  the  vicinity  of  128,  as  is  anticipated,  will 
not  adversely  effect  results. 

There  remains  a  substantial  amount  of  work  to  be  done  in 
assessing  the  performance  of  the  MEM  in  case  #F  >  2.  In  the 
two- frequency  signal  case,  operation  counting  shows  that  Gauss 


elimination  Is  decidedly  superior  to  use  of  Cramer's  Rule. 
Also,  It  Is  likely  that  use  of  some  iterative  algorithm  for 
finding  polynomial  roots  is  better  than  use  of  the  quartic 
formula  (which  generally  requires  use  of  the  cubic  formula) 
for  solving  (2.3)  with  p  =  4.  Of  course,  for  #F  >  3  (p  >  6) 
no  formula  exists  for  solving  (2.3).  A  substantial  portion 
of  the  computing  time  is  required  by  the  COV  subroutine, 
which  is  probably  as  efficient  as  is  possible.  With  the  real¬ 
time  goal  in  mind,  it  is  vital  that  the  most  efficient  algo¬ 
rithms  for  solving  linear  systems  and  polynomial  equations  be 
found . 
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APPENDIX 


CMETA  PROGRAM  META 

DOUBLE  PRECISION  TSA , TPI , TC , Wl 
COMMON/  SAMP/ S  (  20  48) 

COMMON/ FREQ/F  (2048) 

COMMON/RMACJ/R  (  2048) 

COMMON/  PHI/  P  (2,  2) 

COMMON/  PHG/  PO  (2) 

COMMCN/TFRQ/Fl 

C  ASSIGN  VALUES  TO  VARIABLES 

DATA  TSA, TPI/ 8 3. 3 33 333 3333 3333 33D-6, 6. 28318530717958647/ 

DATA  PI,  FC,F  1/  3. 141 59  265, 1909. 8  59  31, 3250./ 

DATA  KU, A, UP, N, AM, SD/ 2 048, 1.5, 2, 3,0. ,0./ 

NN=NP+N 

KL=NN+1 

C  LOOP  TO  CALCULATE  AND  PRINT  PREDICTED  FREQUENCIES  FROM  (4.4)  AND 

C  (4.5) AND  MEAN  AND  STANDARD  DEVIATION  THEREOF  FOR  NUMBER  OF  BITS  = 

C  2, 3,  ...  ,16 

DO  5  NB=2,16 

WRITE  (6,200)  NB 

200  FORMAT  (10H  NUMBER  OF  BITS  =  , 1 3/ / ) 

C  CALCULATE  EXACT  TRANSMITTED  SIGNAL  VALUE  FROM  (4.1) 

Wi  =F  1*  T  PI 
DO  10  J=  1,  2  048 
FJ= J-l 
TC=F  J*TSA 
TS=Wl*TC 

10  S ( J)=A* SIN (TS) 

C  SIMULATION  OF  ADC  QUANTIZATION  OF  THE  SIGNAL 

NBP=  2*  *NB 
PNB=NBP 
DV= 10./ PNB 
DVI=PNB/10. 

DO  15  J=  1,  2  048 
X=S (J) 

X=  X*  D  VI 

KP=X 

X=KP 

15  S(J)  =  X*DV 

CALCULATE  COEFFICIENTS  FOR  AND  SOLVE  LINEAR  SYSTEM  (3.2)  USING 
(4.6)-<4.8) 

DO  20  K=KL,KU 
KK  =  K 

CALL  COV (NP, NN, KK) 

DEL=P  (1,1  )*P  (2,2)-P(l,2)*P  (2,1) 

Al=(P(l,2)*PO(2)-PO(l)*P  (2,2)  )/DEL 
A2=(P(2,l)*PO(l)”PO(2)*P(l,l)  )/  DEL 
C  CALCULATE  AND  PRINT  PREDICTED  FREQUENCY 

F ( K ) =F  C*  A  RCTA ("Al, SQRT (4*a2"A1* Al ) ) 

R(K)=SQRT(A2) 

20  WRITE  (6,  300)  K,R(K),F(K) 

30  .0  FORMAT  (  IX  .  1 8 , 2G20. 8  ) 
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C  CALCULATE  AND  PRINT  MEAN  AND  STANDARD  DEVIATION  OF  PREDICTIONS 

CALL  STATS  ( AM , SD , KL , KU ) 

5  WRITE  (0,400)  AM ,  SD 
400  FORMAT  ( IX, 2C20. b/// ) 

STOP 

END 

SUBROUTINE  C OV (N P , NN , LP ) 

COMMON/ SAMP/ S  (2040) 

COMMON/  PHI/  P  (2,  2  ) 

COMMON/ P  HQ/  PO  (  2  ) 

C  CALCULATE  THE  PHI(I,J)  OF  (3.3) 

C  CALCULATE  DIAGONAL  ELEMENTS,  PHl(J,J),  OF  COVARIANCE  MATRIX  - 

C  ASSIGN  TO  P(J,J) 

L=LP-1 
N  I*  N  N~  N  P 
NL=LP~NI 
B=0. 

DO  5  J=NL,L 
5  B  =  B+S  (J)*S  (J) 

DO  10  J= 1 , NP 
K=LP-J 
I  =  NL“J 

B*B  +  S(I)*S  (I  )-S  (K)*S  (K) 

10  P(J,J)=B 

C  CALCULATE  REMAINDINC  PHI(I,J) 

DO  15  KK= 1 , NP 
B=  0. 

C  CALCULATE  PHI(0,KK)  -  ASSIGN  TO  PO(KK) 

DO  20  J=1,NI 
N=LP-J 
M=N" KK 

20  B=B  +  S  (N)*S  (M) 

PO  (KK)  =  B 

C  CALCULATE  P HI  ( I , J ) , J" 1= KK , 1< I< N P- 1 , 2< J< N P  - 

IF  (KK.EQ.NP)  GO  TO  15 
DO  25  K=1,NP-KK 
I=K 

J=KK+K 
N=LP-K 
M=N“ KK 
Nl=NL-K 
Ml =N 1“ KK 

B=B  +  S (Nl)*S (Ml )“S  (N)*S  (M) 

P  (I  ,  J)=B 
25  P(J,I)=E 

THE  PREVIOUS  STATEMENT  TAKES  ADVANTAGE  OF  SYMMETRY 
COVARIANCE  MATRIX 
15  CONTINUE 
RETURN 
END 


ASSIGN  TO  P(I,J) 


OF 
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CM ETA  PROGRAM  META 

DOUBLE  PRECISION  TSA.TPI , TC,Wl 
COMMON/ SAMP/ S  (20  48) 

COMMON/  F  REQ/  F  (  2  0  48  ) 

COMMON/RMACV  R  (  2048) 

COMMON/  PKI/  P  (2,  2) 

COMMON/  PHQ/  PO  (2) 

COMMON/  TFRU/Fl 

C  ASSIGN  VALUES  TO  VARIABLES 

DATA  TSA.TPI/ 8  3.  3 33 333 33 33 3 333  33d- 6, 8. 2  8  318  53  07179  5864  7/ 

DATA  PI,  FC.Fl/  3. 141 59  265, 1909. 8  59  31, 3250./ 

DATA  KU  ,A,NP,  N,AM,  SIV  2048, 1. 5, 2, 3, 0  .  ,0./ 

NN=NP+N 

KL=NN+1 

C  LOOP  TO  CALCULATE  AND  PRINT  PREDICTED  FREQUENCIES  FROM  (4.4)  AND 

C  (4. 5) AND  MEAN  AND  STANDARD  DEVIATION  THEREOF  FOR  NUMBER  OF  BITS  = 

C  2, 3,  ...  ,16 

DO  5  NB=2,16 

WRITE  (6,200)  NB 

200  FORMAT  (18H  NUMBER  OF  BITS  =  , 1 3/ / ) 

C  CALCULATE  EXACT  TRANSMITTED  SIGNAL  VALUE  FROM  (4.1) 

Wl=Fl*TPI 
DO  10  J=  1 , 2  0  48 
FJ=J-1 
TC=F  J*TSA 
TS=Wl*TC 

10  S ( J)=A* SIN (TS) 

C  SIMULATION  OF  ADC  QUANTIZATION  OF  THE  SIGNAL 

NBP=2**NB 
PNB=NBP 
DV=10./ PNB 
DVI=PNB/  10. 

DO  15  J=  1 ,  2  0  48 
X=S (J) 

X=  X*  D  VI 

KP-X 

X=KP 

15  S(J)  =  X*DV 

CALCULATE  COEFFICIENTS  FOR  AND  SOLVE  LINEAR  SYSTEM  (3.2)  USING 
(4.6)-  (4.8) 

DO  20  K=KL,KU 
KK  =  K 

CALL  COV (NP, NN.KK) 

DEL=P  (1,1)*P  (2,2)-P(l,2)*P  (2,1) 

Al=(P(l,2)*PO(2)-PO(l)*P  (2,2)  )/DEL 
A2=  (P  (2,1  )*PO(l)”PO(2)*  P  (1,1)  )/DEL 
C  CALCULATE  AND  PRINT  PREDICTED  FREQUENCY 

F  (K)=FC*ARCTA  ("Al,  SQRT  (4*A2-Al*Al )  ) 

R(K)=SQRT(A2) 

20  WRITE  (  6,  300)  K,R(K),F(K) 

30.0  FORMAT  (  IX  .  1 8 . 2G  20. 8  ) 
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C  CALCULATE  AND  PRINT  MEAN  AND  STANDARD  DEVIATION  OF  PREDICTIONS 

CALL  STATS  ( AM , SD , K L , KU ) 

5  WRITE  (0,400)  AM,SD 
400  FORMAT  (  IX,  2C20. 6/// ) 

STOP 

END 

SUBROUTINE  C OV (N P , NN , LP ) 

COMMON/ SAMP/ S  (204o) 

COMMON/  PHI/  P  (2, 2  ) 

COMMON/  P  HQ/  PO  (2) 

C  CALCULATE  THE  PHl(I,J)  OF  (3.3) 

C  CALCULATE  DIAGONAL  ELEMENTS,  PHl(J,J),  OF  COVARIANCE  MATRIX  - 

C  ASSIGN  TO  P(J,J) 

L=LP"1 

NI=NN“NP 

NL=LP-NI 

B=0. 

DO  5  J=NL,L 
5  B=B+S  (  J)*S  (J) 

DO  10  J=1,NP 
K=LP”J 
I  =  NL-J 

B=B  +  S  (I)*S(I)-S(K)*S(K) 

10  P(J,J)=B 

C  CALCULATE  REMAINDINC  PHI(I,J) 

DO  15  KK= 1 , NP 
B=  0. 

C  CALCULATE  PHI(0,KK)  -  ASSIGN  TO  PO(KK) 

DO  20  J=1,NI 
N=LP~J 
M=N“ KK 

20  B=B  +  S  (N)*S  (M) 

PO  (KK)  =  B 

C  CALCULATE  P  J!I  ( I  ,  J  )  ,  J"I=  KK,  1<  I<  N  P- 1 ,  2<  J<  N  P  -  ASSIGN  TO  P(I,J) 

IF  (KK.EQ.NP)  GO  TO  15 
DO  25  K-l,  NP-KK 
I  =  K 

J=KK+K 
N=LP~K 
M=N“ KK 
Nl^NL-K 
Ml=N 1-KK 

B=B+S  (Nl)*S  (Ml  )~S  (N)*S  (M) 

P  (I  ,  J)=B 
25  P(J,I)=R 

THE  PREVIOUS  STATEMENT  TAKES  ADVANTAGE  OF  SYMMETRY  OF 
COVARIANCE  MATRIX 
15  CONTINUE 
RETURN 
END 
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SUBROUTINE  STATS  ( AM , S D , KL , KU ) 

STATS  CALCULATES  THE  MEAN  AND  STANDARD  DEVIATION  OF  THE 
KU-(NP+N)  PREDICTED  FREQUENCIES  F(J) 

COMMON/  F  REQ/  F  (204b) 

COMMON/  TFRQ/Fl 

51  =  0. 

52  =  0. 

XN=KU_  KL+  1 
DO  5  J=KL,KU 
Sl=Sl+F  (J) 

5  S  2 “S  2  +  (F  ( J  )  “F 1  )*  *  2 
AM=Sl/XN 
SD=SQRT  (S2/XN) 

RETURN 

END 

FUNCTION  A RCTA ( X, V ) 

ARCTA  CALCULATES  RADIAN  FREQUENCY  DETERMINED  BY  ROOT  X+JY  OF 
(2.3)  BY  USE  OF  (4.4) 

DATA  PI, HP  1/  3. 14159  26  5, 1.5707963  2/ 

IF  (X)  1,2,3 

1  ARCTA=ATAN(Y/X)+PI 
RETURN 

2  ARCTA=  H PI 
RETURN 

3  ARCTA=ATAN  (Y/X) 

RETURN 

END 
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MAXIMUM  ENTROPY  SPECTRAL 


DEMODULATOR  INVESTIGATION  II 

1.  INTRODUCTION 

A  transmitted  signal,  which  is  masked  by  noise  and,  possibly,  interfer¬ 
ence,  is  assumed  to  be  a  sinusoid  with  Hertz-frequency  which  varies  over  some 
finite  set  F  of  positive  real  numbers.  Let  s(t)  be  the  value  of  this  continu¬ 
ous-time  signal  at  time  t.  Of  particular  interest  is  the  case  where  t?F  =  2 ; 
that  is,  where  the  frequency  is  switched  back  and  forth  between  two  distinct 
values . 

The  object  of  this  study  is  to  evaluate  the  performance  of  the  Maximum 
Entropy  Method  (=  MEM)  of  spectral  estimation  for  short  segments  of  the  dis¬ 
crete-time  signal  which  results  from  sampling  s(t)  every  T  seconds.  To  be 
more  specific,  we  wish  to  determine  the  accuracy  of  the  MEM  estimates  of  the 
transmitted  signal  frequencies  for  the  cases  (a)  received  signal  =  trans¬ 
mitted  signal  +  noise  and  (b)  received  signal  =  transmitted  signal  +  noise  + 
interference. 

The  ultimate  goal  of  this  project  is  the  real-time  use  of  the  MEM  (if 
feasible)  to  identify  the  frequencies  of  a  transmitted  signal,  thereby 
countering  the  effects  of  noise  and  interference  (such  as  that  produced  by 
a  keyed,  slewing,  or  CW  jammer). 

In  this  paper,  we  obtain  a  result  (Corollary  6.5)  which  explains  the 
excellent  performance  of  the  MEM  in  a  simulation  [4]  for  the  simple  case  of 
a  single-frequency  sinusoidal  signal  with  no  noise  (s(kT)  =  Asin(2nfkT  +  P)). 
An  effort  to  resolve  the  problem  in  case  (a)  above  (t(kT)  =  s(kT)  +  n(kT) , 
where  s(kT)  is  as  above  and  n(kT)  is  independent,  zero-mean,  Gaussian  noise) 


has  not  been  successful.  This  effort  will  be  continued  with  the  support  of 
Grant  #AF0SR  78  -  3614. 

2 .  THE  MODEL 

We  model  the  sampled  signal  s(kT)  as  the  output  of  a  linear  discrete 
system  (=  LDS)  [1];  in  particular,  by  the  difference  equation 

P 

(2.1)  s(kT)  =  -  l  a.  s((k-j)T)  +  u(kT), 

3=1  : 

where  the  are  real  numbers,  T  is  the  sampling  period,  and  u(kT)  is  the 
transmitted  signal  at  time  kT.  As  (2.1)  enables  one  to  "predict"  s(mT)  from 
s((m-l)T),  ...  ,  s((m-p)T),  and  u(mT) ,  the  name  "linear  prediction"  is  also 
associated  with  the  MEM. 

There  are  several  other  equivalent  representations  of  the  LDS  given  by 
(2.1)  [1;  pp.  85-86].  For  frequency-domain  considerations,  the  representa¬ 
tion 

S(z)  =  H(z)  U(z) , 

where  S(z)  and  U(z)  are  the  z-transforms  of  s(kT)  and  u(kT)  respectively  and 
H(z)  is  a  rational  function  of  z  called  the  "system  transfer  function,"  is 
useful  [1;  pp.  220-  282].  The  transfer  function  corresponding  to  (2.1)  is 
given  by 

P 

H(z)  =  1/  (1  +  a^  z--1). 

Thus  H  is  an  "all-pole"  transfer  function  with  p  poles,  namely,  the  p  solu- 
tions  of  the  equation  1  +  a_.  z~^  =  0  or,  equivalently  (if  z  i  0,  as  is 

required  by  the  definition  of  the  z- transform) , 


3.  INTIMATION  OK  MOD1X,  PARAMITHRS 


'lilt1  nodal  parameters  Uj , . y,  in  (.'.1)  aiv  appixsximated  by  the  usual 

type  ot  least-squuvs  am  lysis  in  the  t  ime-donuin  ( ? ;  pp.  563  -  567] .  Tor  a 
given  positive  integer  m,  we  piwliet  s(n\T)  to  lv 


(3.1) 


»(m'rt  dl  -  y  a.  s((m-  j)T), 
3*1  1 


where  the  a.  are  chosen  go  as  tv'  minimize  tie  mean  square  error 
m-1 

)'  e‘(KT)  (e(kT)  <*f  s(kT)  -  sTkTT) 
k=m-N 


corivsponding  to  the  pivceding  N  s.unples  (theivhv  using  the  "covariance 
method"  ot  linear  prediction  p.  So‘t|).  ITiis  leads  to  the  system  of 
lineal'  equations 

P 

l  a.  (m,  N)  A.  .(m,  N)  =  ,.(m,  N)  (i  =  l,2,...,p), 

j-1  J  ij 


(3.:) 


wliere ,  for  all  (i,  j)  t  (0 . p)  x  { 1 ,  ...  ,p) 


(3.3) 


N-l 

fr(m,N)  d£  l  s«m-N+k-i)T)  s((m-N  +  k-  i)T), 
J  k=0 


toi'  determining  the  a.  which  yield  the  prvdiction  of  s(mT).  Ptom  (3.3)  and 
(3.3),  we  see  that.  N  ♦  p  consecutive  samples  s((m-N-p)T),  ...  ,  s((m-  1)T) 

'Uv  needed.  ‘lints  if  we  do  tot  allow  negative  arguments ,  s((N  +  p+l)T)  is  the 
first  sample  we  con  piwlict  . 

In  making,  the  predict  ion  (3.1),  we  are  assuming’,  that  u(mT)  is  eom- 
pletely  unknown,  which  is  often  the  case. 


11-5 


a 


a. _ ' -  . . . 


■MMtei 


4 .  TOE  ANALYTICAL  APPROACH 


The  determination  of  the  frequencies  f^  e  F  exhibited  by  the  trans¬ 
mitted  signal  involves  (a)  calculating  the  <J>Mm,N)  from  (3.3),  (b)  solving 
the  system  of  linear  equations  (3.2)  for  the  a.  ,  (c)  solving  the  polynomial 
equation  (2.2)  with  "a/'  in  place  of  "a/',  (d)  determining  the  radian- 
frequency  0^  in  the  polar  form  r^  exp ( ± i0 . )  of  each  of  the  complex  conjugate 
pairs  of  roots  of  (2.2),  and  (e)  multiplying  6 ^  by  an  appropriate  real 
number  to  obtain  the  Hertz- frequency. 

The  problem  of  determining  a  confidence  interval  for  f^ ,  which  requires 
finding  a  probability  density  function,  is  difficult,  even  for  the  simple 
case  F  =  {f^}  with  noise.  This  case,  in  which  a  two-pole  model  (p  =  2  and 
the  equation  (2.2)  is  quadratic)  is  appropriate,  was  considered  by  the 
writer  in  1977  under  the  USAF/ASEE  Sunmer  Faculty  Research  Program.  Late  in 
this  program,  this  analytical  effort  was  abandoned  in  favor  of  a  computer 
simulation.  In  this  paper,  we  return  to  the  theoretical  effort. 


5.  COMPUTER  SIMULATION 

In  the  above-cited  simulation,  we  assumed  F  =  {f^}  and  the  transmitted 
signal  is  A sin(2irf^T) ,  with  f^  e  (0,  6000)  and  T  =  1/12000  (in  accordance 
with  The  Sampling  Theorem  [1;  p.  291]),  ignoring  noise  of  the  type  cited  in 
[3;  pp.  3-5]  and  interference;  thus  the  only  noise  is  "quantization  noise" 
(from  use  of  a  simulated  analog- to-digital  converter  (=  ADC)  in  sampling  the 
continuous -time  signal)  and  "roundoff  noise"  (from  performing  the  MEM  calcu¬ 
lations  with  a  digital  computer) .  This  simulation  shewed  excellent  perform¬ 
ance  of  the  MEM  for  sampling  with  12  -  16  bit  ADCs  (which  are  available 
commercially)  if  f^  is  not  too  close  to  either  0  or  6000  [4;  pp.  10-26]. 


The  table  below  from  (4;  p.  17]  gives  the  arithmetic  mean  (=  AM)  and 
standard  deviation  (=  SD)  of  2043  predicted  frequencies  for  the  number  of 
error  terns  used  in  the  minimization  N  =  2,  3,  ...  ,12  with  A  =  1.5,  f  =  3250, 
and  simulation  of  a  16-bit  ADC. 


N 

AM 

SD 

2 

3249.9999 

0.045065880 

3 

3250.0008 

0.034125918 

4 

3249.9995 

0.024391433 

5 

3249.9991 

0.021149200 

6 

3250.0010 

0.019371934 

7 

3249.9993 

0.015833400 

8 

3249.9985 

0.013296252 

9 

3249.9994 

0.011714947 

10 

3249.9999 

0.010651756 

11 

3249.9995 

0.0083492643 

12 

3250.0000 

0.0071696392 

6.  THE  TWO-POLE  CASE  WITHOUT  NOISE 

In  this  section,  we  give  some  theoretical  results  which  explain  the 
excellent  performance  of  the  MEM  in  the  above-cited  simulation 

Throughout  the  remainder  of  this  paper,  we  use  R  and  Z+  to  denote  the 
set  of  real  numbers  and  the  set  of  positive  integers  respectively. 

6.1.  Lenina.  For  all  B,  P  e  R,  k,  m,  n  e  Z+, 

sin((k  +  m)B  + P)  sin((k  +  n)B  +  P)  -  sin((k  +  m  +  n)B  +  P)  sin(kB  +  P) 

=  sin(mB)  sin(nB) . 

Proof.  By  familiar  trigonometric  identities,  the  left-hand  side  of  the 
above  equality  is  equal  to 

(1/2)  (cos((m- n)B)  -  cos((2k  +  m  +  n)B  +  2P) ) 

-  (1/2)  (cos((m  +  n)B)  -  cos((7k  +  m  +  n)B  +  2P)) 
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=  (1/2)  (cos((m-n)B)  -  cos((m  +  n)B) 


=  sin(  ((m  +  n  +  m  -  n)/2)B)  sin(( (m  +  n  -  (m  -  n))/2)B) 


sin(mB)  sin(nB) . 


6.2.  Corollary.  If,  for  all  k  e  Z  , 


<•%.  s(kT)  =  A  sin(2irfkT  +  P) , 

where  f ,  P,  T  e  R,  then,  for  all  k,m,  n  e  Z+ 

s,  ,  s,  ,  -  s,  ,  ,  s.  =  A2  sin(2irfmT)  sin(2?rfnT) . 

k+m  k+n  K+m+n  k 

In  particular,  for  m  =  n  =  1,  we  have  +  ^2  -  +  2  =  A2  sin2(2TrfT). 

6.3.  Theorem.  Ifp  =  2,  N  e  Z+-{1},  m  =  N  +  p  +  1,  is  defined  as  in 
Corollary  6.2,  0  <  P  <,  2it,  and  0  <,  f  <_  1/(2T),  then  (a)  the  system  (3.2)  is 
singular  if  and  only  if  f  =  0  or  f  =  1/(2T)  and  (b)  in  the  non-singular  case 
its  solution  is  the  pair  (a^,  a^)  given  by  a^"  =  -2cos(2irfT)  and  aj  =  1. 

Proof,  (a)  If  p  =  2  and  m  =  N  +  p  +  1,  then  (from  (3.3)) 


.j(m»  N)  =  jJ0  s3  +  k-  is3  +  k-  j 


and  the  system  (3.2)  is  as  follows: 

(s22  +  s32  +  ...  +  sN+12)  ax  +  (Bj  s2  +  s2  s3+...+sN  sN+1)a2 


(6.1) 


=  -  (s2  s3  +  s3  s4  +  . . .  +  sN+1  sN+2) 
(S;LS2  +  S2  s3+  ...  +  sNsN+1)a^  +  (s-j2  +  s22  +  ...  +sN2)i^ 


=  -  (s,  s„  +  s„  s„  +  . . .  +  s,,  s„.  „)  . 


Let  D  be  the  determinant  of  the  matrix  of  coefficients  in  (6.1).  Then 

D=  <Sl2*‘"t3N2)(s22*---*W)  -  <slVs2s3*---tsNW: 

It  is  convenient  to  arrange  the  terms  of  D  in  a  2N  by  N  array 


where  both  U  and  V  are  N  by  N  arrays ,  u . .  df  s .  2  s .  2  and  v  df 

x]  i  j  + 1  i  j  — = 

Si  si  +  1  sj  s  j  +  1  ^01  ajT  ^ »  i  ^  e  {1,2,...  ,N}~.  The  main  diagonal  ele¬ 
ments  of  U  cancel  with  the  corresponding  elements  of  V;  that  is,  u  +  v 

ii  ii 

=  0  for  all  i  e  {1,2,  ...  ,N}.  The  remaining  terms  of  D  can  be  grouped  as 
follows : 

<Tij  “ }  Uij  +  Vij  +  uji  +  vji  (i  =  1,  . . .  ,N-1;  j  =  i  +  l,  ...  ,N). 


N-l  N 

D  =  l  l  T. .. 
i=l  j=i+l 


It  is  helpful  to  replace  this  sum  by  one  in  which  we  sum  along  "lines" 
j  ~  i  +  k »  k  =  1,2,...  ,  N -  1 ;  that  is 

N-l  N-k 


D  =  I  It.  . 

k=l  i=l  1  k 


Obviously , 


Tij  =  Si2sj+12  "  sisi+lsjsj+i  +  sj2si+i2  ~  sjsj+isisi+l 
=  (sj  si  +  l  ~  sj  +isi)2* 

Hence,  by  use  of  Corollary  6.2  (with  k  =  i,  m  =  k,  and  n  =  1) ,  we  have 
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(A2  sin(2irfkT)  sin(2TrfT))2 


A4  sin2(27rfT)  sin2(2irfkT) 


D=  l  (N-k)  A4  sin2(2irfT)  sin2(2irfkT) 
k=l 


A4  sin2(2TrfT)  £  (N  -  k)  sin2(2irfkT) 


Hence,  D  =  0  if  and  only  if  sin(2irfT)  =  0  or  sin(2irfkT)  =  0  for  all  k  e 
{1,  2,  ...  ,  N-  1},  which  is  equivalent  to  sin(2irfT)  =0.  As  0  <,  f  <  1/(2T) , 

0  <  2irfT  <  *  and,  therefore,  sin(2TrfT)  =  0  if  and  only  if  f  =  0  or  f  =  1/(2T) 


(b)  In  case  D  i  0,  Cramer's  rule  gives  the  solution  (I^,a2),  where 

^(1)  /r,  —  _  ~(2)  ,  ^ 


^+1^+2 )  (sis2  +  s2  s3  + 


•  ■  . . . . .  ■  — »— 


I 


hi™  «t“il<1>**«(U-Ji(1,*vli<1>  Wl... 

It  is  easy  to  show  that 

Tij(1>  =  lSiSj*l  -  ajW(si*2sj  -  sisj*2>- 

New,  by  two  uses  of  Corollary  6.2,  we  have 

Ti,  i  +  k  )  =  (si  si  +  k  +  1  "  si  +  k  si  +  1}  (si  +  2  si  +  k  si  si  +  k  +  2) 
=  (-  A2  sin(2TtfkT)  sin(2nfT))  (A2  sin(2irfkT)  sin(2nf2T)) 


A4  sin(2irfT)  sin(4irfT)  sin2(?irfkT) . 


Thus 


N-l  N-k 


SD 


y  y  t.  . 

k=l  i=l  1  K 


(1) 


N-l 


=  -  A4  sin(2irfT)  sin(4irfT)  £  (N  -  k)  sin2(2rrfkT) . 

k=l 


(2) 


riially,  let  u. . 21  -  ^  ^  Sj  Sj  ,  2,  Vy*2’  df  s.  +  l  +  ,  s.  s.  + 


i  j  +  l’ 


and 


T..<2)  df  u..(2)  +v..(2)  +u..(2)  +  v..(2) 

ID  —  ID  ID  Di  D1 

(si  +  1  sj  '  si  sj  +  l5  (si  +  2  sj  +  1  "  si  +  1  sj  + 


By  Corollary  6.2, 


(2) 


Ti,  i  +  k  "  (si  +  1  si  +  k  "  si  si  +  k  +  l5  (si  +  2  si  +  k  +  1  "  si  +  1  si  +  k  +  2J 
(si+lsi+k  "  sisi+k+l)  (s(i+l)+l  s(i+l)+k  "  si+l  s(i+l)+k+l) 
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.  ;■  -c.;-  C ,  '  V 


-  '  -  :^Au  dfcAiaU  gV  ,  • 


. mi  III  mill  I  nil  I  M 


=  (A2  sin(2nfT)  sin(2nfkT))  (A2  sin(2irfT)  sin(2ufkT)) 


=  A4  sin2(27rfT)  sin2(27tfkT) 
=  Ti,  i  +  k* 


(2) 

Hence,  Dv  ' 
Thus, 


D. 


D 


-A4  sin(2irfT)  sin(4rrfT)  (N-k)  sin2(2irfkT) 
A4  sin2(2irfT)  (N-k)  sin2(2TrfkT) 


and 


=  -  2  cos(2irfT) 


We  should  remark  that  D  could  have  been  evaluated  easily  by  use  of 
Lagrange's  identity,  however,  this  is  not  the  case  for  and  D(2\ 

6.4.  Remarks.  In  case  p  =  2,  the  minimization  of  Section  3  must 
involve  at  least  two  error  terms;  hence,  we  have  assumed  N  i  1  in  Theorem  6.3. 
In  case  N  =  1,  the  system  (6.1)  takes  the  form 


S1  s2  a2 


S2  S3 


S1  S2  al  + 


a2  =  "  S1  s3 


and  D  =  =  0.  Thus,  (6.1)  does  not  have  a  unique  solution. 

There  is  no  loss  of  generality  in  assuming  m  =  N  +  p+1  in  Theorem  6.3. 
This  is  an  obvious  consequence  of  Corollary  6.2.  For  example ,  ifm  =  N  +  p+  2, 
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then  each  of  the  subscripts  in  the  coefficients  of  the  system  (3.2)  and  in 
the  expressions  for  D,  and  ‘  ^  is  one  larger  than  in  the  case  of 

m  =  N  +  p  +  1;  thus,  by  Corollary  6.2,  D,  and  have  the  same  values. 

Examination  of  the  cases  N  =  2  and  N  =  3  suggested  the  proof  of 
Theorem  6.3.  It  may  be  helpful  to  give  the  proof  of  this  theorem  in  case 
N  =  2.  First, 


D  =  (s^2  +  s?2)  (s?2  +  s32) 


(sls2  +  s2s3)? 


S1  S2S1  s2 


S1  S2  S2  s3 


-  s2  s3  s.^  s2 


s2  S3  s2  s3 


C  2  q  2  |  q  4 

S1  S3  S2 


S1  S22  S3  ~  S1  S22  S3 


=  (s22  -  s1s3)2 

=  (A2  sin2(2irfT))2 
=  A4  sin4(2nfT) . 

Obviously,  D  =  0  if  and  only  if  sin(2nfT)  =  0  or  if  and  only  if  f  =  0  or 
f  =  1/C2T).  Also, 

=  -  (S^2  +  S02)  (S.,  Sj  +  S3  S4)  +  (s^  S3  +  Sp  S^)  (s^  S,  +  S,  S.?) 
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-  S12s2S3  "  sl2s3s4 

-  s22  s2  s3  -  s2-  s3  s4 

+  S1S3S1S2  +  S1S3S2S3 
+  s2  s4  sx  s2  +  s2  s4  s2  s3 

'  V  s3  s4  ~  S23  S3  +  S1  S2  S32  +  S1  S22  S4 

-  (s22  -  Sg)  (s2  S3  -  S4) 

-  (A2  sin2(2irfT))  (A2  sin<2irfT)  sin(2irf2T)) 
-A4  sin3(25rfT)  sin(4fffT) 


-  (so2  +  S3;  )  (sx  S3  +  S4)  +  (s3  S2  +  So  Sg)  (s2  S3  +  s3  s4 ) 

V  s1  s3  -  s22  s2  s4 

-  s32  S1  s3  -  s32  s2  s4 

+  S1  S2  s2  S3  *  S1S2S3S4 
+  s2  s3  s2  s3  +  s2  s3  s3  s4 

-s23s4  -  Sls33  +  s1s2s3s4  +  s22s32 

(s22  -  9ls3)  (s32  -  s2s4) 
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=  (A2  sin2(2*fT))  (A2  sin2(2wfT)) 


=  A4  sin4(2irfT) 

=  D. 

Hence,  if  D  i  0, 

_  D(1)  -A4  sin3(27rfT)  sin(4rrfT) 

a,  =  -  =  -  =  -2cos(27rfT) 

D  A4  sin4(27rfT) 


and 


D 

D 


1. 


=  -  (- 2cos(2irfT))  (A  sin(2irf  (N  +  p)T  +  P) ) 

-  (1)  (A  sin(2itf  (N  +  p  -  1)T  +  P) ) 

=  A(sin(2irfT  +  2irf(.N  +  p)T  +  F>  -  sin(277fT- 27rf(N  +  p)T  -  P)) 

-  A  sin(2?rf  (N  +  p  -  1)T  +  P) 

=  A  sin(  2irf  (N  +  p  +  1)T  +  P)  +  A  sin(  2irf  (N  +  p  -  1)T  +  P) 

-  A  sin(2irf  (N  +  p  -  1)T  +  P) 

=  A  sin(2itf  (N  +  p  +  1)T  +  P)  =  s((N  +  p  +  l)T). 


6.6.  Remark.  Corollary  6.5  explains  the  remarkably  good  simulation 
results  cited  in  Section  5  for  the  case  of  no  noise.  We  assumed  in  Theorem 
6 . 3  and  Corollary  6 . 5  that  the  f irst  N  +  p  (=  N+2)  samples  of  the  received 
signal  are  actually  the  first  N+p  samples  of  the  transmitted  signal  and 
shewed  that  the  minimization  process  of  linear  prediction  by  the  covariance 
method  leads  to  values  of  a-^  and  a2  such  that  the  predicted  signal  sample 
s((N  +  p  +  1)T)  is  equal  to  the  transmitted  signal  sample  s((N  +  p+l)T).  Thus 
there  is  no  model  error,  according  to  this  corollary;  all  the  error  is  due 
to  the  quantization  and  the  computational  process. 

6.7.  Theorem.  Suppose  the  hypotheses  of  Theorem  6 . 3  hold  with 
0  <  f  <  1/(2T)  (so  that  the  system  (3.2)  is  non-singular). 

(a)  The  zeros  of  z2  +  z  +  a2  (that  is,  the  poles  of  H(z)),  with 
a^  =  -  2  cos(2TrfT)  and  =  1  as  in  Theorem  3.2  (b) ,  are 

z1  df  cos(2irfT)  -  .i  sin(2nfT)  =  exp(-  2irfTi_) 


z2  =  cos( 2irfT)  +  jLsin(2iTfT)  =  exp(2*fTi) . 

(b)  These  zeros,  z^  and  z? ,  have  magnitude  r  =  1  (that  is,  can  be 
represented  by  points  on  the  unit  circle)  or,  in  terms  of  a^  and  aj, 

r  =  /  a„  . 


(c)  If  T  =  1/12000  (so  that  0  <  f  <  6000),  then 


(6.2)  f  =  ^  3000 


(6000  /  it)  tan  (Im(z0)  /  Re(z0) ) 

L  L 


(if  0  <  27rfT  <  tt/2), 
(if  27rfT  =  ir/2). 


(6000/n)  (tan“-L(Im(z2)  /Re(z2))  +  tt)  (if  tt/2  <  2TrfT  <  it) 
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or,  in  tents  of  and  a,,, 


(6.3)  f  =  < 


(6000/tt)  tan'^C-  v  4a„-a^2/a^) 


3000 


C 6000/it )  (tan-^(-  »  4  a_,  -  a^2  /  a^)  +  it) 


(if  a^  <  0) , 
(if  a^  =  0), 
(if  a^  >  0) . 


Proof,  (a)  From  Theorem  6.3  (b)  and  the  quadratic  formula,  the  zeros 
of  z2  +  a^  z  +  a",  are 

(-1^/2)  ±  (  /l^2  -  4i;  /  2) 
or,  as  a^2  -  4  =  4cos2(2nfT)  -  4  <0, 

(-  /  2)  ±  i  (  /  4  r,  -  a^2  /  2) 

=  cos(2irfT)  ±  i_sin(2iTfT) . 

(b)  Obviously, 

|  z^  |  =  |  z^  j  =  J (-  a^  /  2)2  +  (  / 4  a~,  -  a^2  /  2)2  =  /a^,  =  J\  -  1 . 

(c)  Let  X  df  Re(z ,)  =  cos(2irfT)  and  Y  df  Im(z9)  =  sin(2irfT). 

Case  1  (0  <  2TrfT  <  it/2).  In  this  case,  X  >  0,  Y  >  0,  and  Y/X  = 
tan(2TifT)  >  0;  thus,  0  <  tan”^(Y  /  X)  <  tt/2  and,  hence,  2TTfT  =  tan  ^(Y/X). 
As  T  =  1/12000, 

f  =  (6000  /  it)  tan-1(Y  /  X)  . 

Case  2  (2TrfT  =  n/2) .  As  T  =  1/12000,  f  =  3000  follows  itmvediately 
frcm  the  case  assumption. 


Case  3  (tt/2  <  2irfT  <  n) .  In  this  case ,  X  <  0 ,  Y  >  0 ,  and  Y/X  = 


tan(2nfT)  <  0;  thus,  -  tt/2  <  tan-1(Y/X)  <  (1  and,  hence,  2nfT  =  tan-1(Y/X)  + 
As  T  =  1/12000, 

f  =  (6000  /  *)  (tan-1  (Y  /X)  +  n  ) . 

Finally,  (6.3)  results  iiimediatoly  from  (6.2),  as  Ro(:"„.,)  =  -  aj  /  2  and 
Im(  :•,.,)  =  <^4  a  ,  -  ,Tj;'/2. 

6.8.  Remarks,  'theorem  6.7  shows  that  the  frequency  of  a  sinusoidal 
signal  without  noise  ot  .uiy  kind  (including  quantization  noise  and  roundoff 
noise)  can  be  determined  (hv  use  of  (6.2)  or  (6.3))  from  a  pole  of  the  all- 
pole  tiunsfer  function  associated  with  the  IJTS  (2.1)  with  p  =  7.  In 
subsequent  work,  we  will  use  (6.3)  to  approx invite  f  in  the  case  of  a  sinus¬ 
oidal  signal  with  tViussian  noise. 

The  simulation  results  given  in  Section  5  show  that  the  stain  Lux  I  devia¬ 
tion  ot  the  pix'dicted  values  ol  f  docivases  as  N  increases .  This  might  seem 
strange  in  view  of  the  fact  that  a  ^  (=  -  2  cos( 2irtT))  and  a.,  (=  1)  are  inde¬ 
pendent  ol  N.  However,  the  computer  program  used  calls  for  the  computation 
ot  p,  r^^.  and  1  which  do  depend  on  N;  that  is,  no  cancellation  is  done 
in  the  computer  computation. 
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MAXIMUM  ENTROPY 

SPECTRAL  DEMODULATOR  INVESTIGATION  III 

1.  INTRODUCTION 

A  transmitted  signal,  which  is  masked  by  noise  and,  possibly,  interfer¬ 
ence,  is  assumed  to  be  a  sinusoid  with  Hertz-frequency  which  varies  over  some 
finite  set  F  of  positive  real  numbers.  Let  s(t)  be  the  value  of  this  continu¬ 
ous-time  signal  at  time  t.  Of  particular  interest  is  the  case  where  #F  =  2; 
that  is,  where  the  frequency  is  switched  back  and  forth  between  two  distinct 
values. 

The  object  of  this  study  is  to  evaluate  the  performance  of  the  Maximum 
Entropy  Method  (=  MEM)  of  spectral  estimation  for  short  segments  of  the  dis¬ 
crete-time  signal  which  results  from  sampling  s(t)  every  T  seconds.  To  be 
more  specific,  we  wish  to  determine  the  accuracy  of  the  MEM  estimates  of  the 
transmitted  signal  frequencies  for  the  cases  (a)  received  signal  =  trans¬ 
mitted  signal  +  noise  and  (b)  received  signal  =  transmitted  signal  +  noise  + 
interference. 

The  ultimate  goal  of  this  project  is  the  real-time  use  of  the  MEM  (if 
feasible)  to  identify  the  frequencies  of  a  transmitted  signal,  thereby 
countering  the  effects  of  noise  and  interference  (such  as  that  produced  by  a 
keyed,  slewing,  or  CW  jammer). 

The  problems  stated  in  the  second  paragraph  appear  to  be  theoretically 
intractable.  We  give  some  partial  theoretical  results  and  discouraging 
simulation  results  for  the  case  (a).  For  case  (b)  without  noise,  we  give 
simulation  results  which  point  to  the  need  for  analog-to-digital  conversion 
of  greater  precision  than  is  possible  with  currently  available  equipment. 


2.  THE  MODEL 


We  model  the  sampled  signal  s(kT)  as  the  output  of  a  linear  discrete 
system  (=  LDS)  [1];  in  particular,  by  the  difference  equation 


(2.1) 


P 


s(kT)  =  -  l 

j-1 


aj  s((k-  j)T)  +  u(kT), 


where  the  a^  are  real  numbers,  T  is  the  sampling  period,  and  u(kT)  is  the 
transmitted  signal  at  time  kT.  As  (2.1)  enables  one  to  "predict"  s(mT)  from 
s((m-l)T),  ...  ,  s((m-p)T),  and  u(mt),  the  name  "linear  prediction"  is  also 
associated  with  the  MEM. 

There  are  several  other  equivalent  representations  of  the  LDS  given  by 

(2.1)  (1;  pp.  85-86].  For  frequency-domain  considerations,  the  representa¬ 
tion 

S(z)  *  H(z)  U(z), 

where  S(z)  and  U(z)  are  the  z-transforms  of  s(kT)  and  u(kT)  respectively  and 
H(z)  is  a  rational  function  of  z  called  the  "system  transfer  function,"  is 
useful  (1;  pp.  220  -  282].  The  transfer  function  corresponding  to  (2.1)  is 
given  by 

P 

H(z)  =  1/  (1  ♦  l  a  z"J). 

j  =  l  J 

Thus  H  is  an  "all-pole"  transfer  function  with  p  poles,  namely,  the  p  solu- 

p  * 

tions  of  the  equation  1  ♦  a.  =  0  or,  equivalently  (if  z  4  0,  as  is 

required  by  the  definition  of  the  z-transform) , 

(2.2)  zp  +•  al  zp_1  ♦  ...  ♦  a  j  z  +  ap  =  0. 
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3.  ESTIMATION  OF  MODEL  PARAMETERS 


The  model  parameters  a^,  ...  ,  a^  in  (2.1)  are  approximated  by  the  usual 
type  of  least-squares  analysis  in  the  time-domain  [2;  pp.  563  -  567].  For  a 
given  positive  integer  m,  we  predict  s(mT)  to  be 


_ P 

(3.1)  s(mT)  df  _  l  a.  s((m-j)T), 

j  =  l  J 


where  the  a^  are  chosen  so  as  to  minimize 


ra-1 

l  e2  (kT)  (e(kT)  df  s(kT)  -  F(kT)) 
k*m-N 


corresponding  to  the  preceding  N  samples  (thereby  using  the  "covariance 
method"  of  linear  prediction  [2;  p.  564]).  This  leads  to  the  system  of 
linear  equations 

P  __ 

(3.2)  ^  aj  (m,  N)  (m,  N)  *  -*0.(ra,N)  (i  =  1,  2,  ...  ,  p) , 

where,  for  all  (i, j)  e  (1,  ...  ,  p}2, 

N-l 

C3'3)  4>..(m,N)  df  £  s  ((m  -  N  ♦  k  -  i)T)  s((m  -  N  +  k  -  j)T) , 

J  k=0 

for  determining  the  a^  which  yield  the  prediction  s (mT) .  From  (3.2)  and 

(3.3) ,  we  see  that  N  +  p  consecutive  samples  s((m  -  N  -  p)T) ,  ...  ,s((m-l)T) 
are  needed.  Thus  if  we  do  not  allow  negative  arguments,  s((N  +  p+  1)T)  is 
the  first  sample  we  can  predict. 

In  making  the  prediction  (3.1),  we  are  assuming  that  u(mT)  is  com¬ 
pletely  unknown,  which  is  often  the  case. 
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4.  THE  ANALYTICAL  APPROACH 


The  determination  of  the  frequencies  f.  c  F  exhibited  by  the  transmitted 
signal  involves  (a)  calculating  the  4>„(m,  N)  from  (3.3),  (b)  solving  the 
system  of  linear  equations  (3.2)  for  the  a y  (c)  solving  the  polynomial  equa¬ 
tion  (2.2)  with  "a 7"  in  place  of  "a/',  (d)  determining  the  radian-frequency  0.. 
in  the  polar  form  r^  exp (±i_0^)  of  each  of  the  complex  conjugate  pairs  of 
roots  of  (2.2),  and  (e)  multiplying  0^  by  an  appropriate  real  number  to  obtain 
the  Hertz- frequency. 

The  problem  of  determining  a  confidence  interval  for  the  frequency  of  the 
transmitted  signal  at  a  given  instant  and  for  the  frequency  of  the  interfer¬ 
ence  (if  any),  which  requires  finding  probability  density  functions,  is 
difficult,  even  for  the  simple  case  F  =  {fj}  with  noise.  This  case,  in  which 
a  two-pole  model  (p  =  2  and  the  equation  (2.2)  is  quadratic)  is  appropriate, 
was  considered  by  the  writer  in  1977  in  the  USAF/ASEE  Summer  Faculty  Research 
Program.  The  analytical  approach  was  abandoned  in  favor  of  a  computer  simula¬ 
tion.  In  this  paper  (Section  6),  we  return  to  this  theoretical  effort. 

5.  SIMULATION  IN  THE  TWO-POLE  CASE  WITHOUT  NOISE 

In  the  above-cited  simulation,  we  assumed  F  =  (f^)  and  the  transmitted 
signal  is  A  sin^irfjT) ,  with  fj  e  (0,6000)  and  T  =  1/12000  (in  accordance 
with  The  Sampling  Theorem  [1;  p.  291]),  ignoring  noise  of  the  type  cited  in 
[4;  pp.  3-5]  and  interference;  thus  the  only  noise  is  "quantization  noise" 
(from  use  of  a  simulated  analog-to-digital  converter  (*  ADC)  in  sampling  the 
continuous-time  signal)  and  "roundoff  noise"  (from  performing  the  MEM  calcu¬ 
lations  with  a  digital  computer).  This  simulation  showed  excellent  perform¬ 
ance  of  the  MEM  for  sampling  with  12-  16  bit  ADCs  (which  are  available 
commercially)  if  fj  is  not  too  close  to  either  0  or  6000  (s»  pp.  10*  26]. 
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The  table  below  from  [5;  p.  17]  gives  the  arithmetic  mean  (*  AM)  and 
standard  deviation  (=  SD)  of  2043  predicted  frequencies  for  the  number  of 
error  terms  used  in  the  minimization  N  *  2,  3,  ...  ,12  with  A  =  1.5,  fj  =  3250, 
and  simulation  of  a  16-bit  ADC. 


N 

AM 

SD 

2 

3249.9999 

0.045065880 

3 

3250.0008 

0.034125918 

4 

3249.9995 

0.024391433 

5 

3249.9991 

0.021149200 

6 

3250.0010 

0.019371934 

7 

3249.9993 

0.015833400 

8 

3249.9985 

0.013296252 

9 

3249.9994 

0.011714947 

10 

3249.9999 

0.010651756 

11 

3249.9995 

0.008349264 

12 

3250.0000 

0.007169639 

6.  THE  TWO- POLE  CASE  WITH  NOISE 

In  this  section,  we  consider  the  case  in  which  the  number  of  poles  p  =  2 


and 


(sk  d|)  s(kT)  =  t(kT)  +  n(kT), 


where  t(kT)  =  Asin(2irfT  ♦  P)  and  n(kT)  is  independent,  zero-mean,  Gaussian 
noise. 

We  assume  for  simplicity  that  N  =  2,  in  which  case  the  system  (3.2)  is 
(s22  +  s32)  a^  +  (Sj  s2  +  s2  s3^  *2  =  "  ^s2  S3  +  S3S4^’ 

(6.1) 

(sl  s2  *  s2  s3^  al  +  ^sl2  +  s22)  a2  =  "  (S1  s3  +  s2  s45  * 


As  was  shown  in  [6;  p.  6,  Theorem  6.3],  this  system  of  equations  is  singular 
if  and  only  if  f  »  0  or  f  *  1  /  (2T),  in  case  we  restrict  f  to  [0,  1  /  (2T) ] . 


In  the  non-singular  case,  its  solution  is  the  ordered  pair  (8j ,  a,)  given  by 


_  ‘  (S12  +  S22)  (s2  S3  +  S3  S4}  +  (sl  s3  *  s2  SA)  (si  s2  +  S2  s3> 

(6.2)  *  _ _ 

(sl2  +  s22)(s22  ♦  s32)  -  (Sj  s,  ♦  s2  s3)2 
=  (s22  -  sj  s3)  (Sj  s4  -  s2  s3)  /  (s22  -  Sj  s3)2 
=  (Sj  s4  -  s2s3)/  (s22  -  SjS3) 


•• 


and 


(6.3) 


(s  2  «•  s  2)(s 


1  S3  *  s2  S4> 


(s,s2 


*  S2  s3)(s2  s3  *  s. 


3  4J 


(Si2  +  S22)(S22  +  S32)  '  (S1S2  +  S2S3)2 
=  (s22  -  Sj  s3)  (s32  -  s2  s4)  /  (s22  -  sx  s3)2 

^S3"  "  s2  s4^  ^S22  "  S1S3^‘ 


The  problem  of  determining  the  probability  density  functions  of  the  ran¬ 
dom  variables  a^  and  a2  is  intractable  according  to  my  colleague  T.  S.  Bolis, 
an  expert  in  probability  theory.  The  numerators  and  (common)  denominator  in 
the  expressions  for  a^  and  a.,  have  distributions  which  look  "something  like" 
non-central  chi-square  distributions.  However,  Bolis  has  shown  (by  use  of  the 
characteristic  function  -  a  complex  variable  analog  of  the  moment  generating 
function)  that  the  numerator  and  denominator  are  dependent  for  both  a^  and  a 
thus  we  cannot  easily  calculate  bounds  on  a^  and  a2  (through  bounding  the 
numerators  and  denominator  -  which  requires  independence) . 

As  shown  in  [6;  pp.  14-  16,  Theorem  6.7],  for  the  two-pole  case  without 
noise  (s(kT)  =  Asin(2nfT  +  P)),  the  frequency,  f,  of  the  transmitted  signal 
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v »' 

v  r  s 


-v'1N--c.r  . *  .V 

'  ■  -  j..  As* 


is  given  by 


% 


(6.4) 


f  = 


(6000/n)  tan_1(-  /at  -  r,2/^) 

3000 

(6000/ it)  (tan_1(-  A 

a2  ~  a!2  /  ^3  +  ») 


(if  aj  <  0), 
(if  al  =  0), 
(if  aj  >  0) . 


We  also  use  (6.4)  to  "predict"  f  in  the  two-pole  case  with  noise.  If  we 
could  bound  a^  and  a.,,  we  could  try  to  get  bounds  on  f  as  given  by  (6.4). 

We  can,  however,  obtain  the  expected  value  and  variance  of  the  numerator 
and  denominator  in  the  expressions  for  both  aj  and  a 2  (Theorem  6.4). 

6.1.  Lemma.  For  all  B,  P  e  R,  k,  m,  n  e  Z+, 

sin((k  ♦  m)B  +  P)  sin((k  +n)B  +  P)  -  sin((k  +  m  +  n)B  +  P)  sin(kB  +  P) 

=  sin (mB)  sin(nB) . 

Proof.  By  familiar  trigonometric  identities,  the  left-hand  side  of  the 
above  equality  is  equal  to 

(1/2)  (cos((m-r.)B)  -  cos ((2k  +  m  +  n)B  +  2P) ) 

-  (1/2)  (cos((m  +  n)B)  -  cos ((2k  ♦  m  +  n) B  +  2P)) 

=  (1/2)  (cos((m-n)B)  -  cos((m  +  n)B) 

=  sin(((m  +  n  ♦  m  -  n)/2)B)  sin(  ( (m  +  n  -  (m  -  n))/2)  B) 

=  sin(mB)  sin(nB)  . 

6.2.  Corollary.  If,  for  all  k  e  Z* , 

(sk  df)  s(kT)  =  Asin(2irfkT  +  P), 
where  f,  P,  T  c  R,  then,  for  all  k, m,  n  e  Z+, 


sk  +  m  sk  +  n  "  sk  +  m  +  n  sk  A2  sin(2»fmT)  sin(2nfnT). 

In  particular,  for  m  =  n  =  1,  we  have  sk  +  ^  -  sk  sk  +  9  =  A2  sin2  (2irfT)  . 

6.3.  Lemma.  If  X  is  a  normally  distributed  random  variable,  E(X)  =  0, 
and  Var(X)  =  a2,  then,  for  all  n  e  Z+, 

'Jr*  +  1 

(a)  E(X"n  )  =  0 
and 

n 

(b)  E(X2n)  =  o2n  TT  (2i  -  1). 

i=l 

Fi-oof.  This  is  a  well-known  result;  (a)  follows  immediately  from  the 
fact  that  if  f  is  an  odd  function,  then  /  f(x)  dx  =  0  and  (b)  follows 
readily  from 

/+°°  x2ne'ax2  dx  =  (2a) _n  (w/a)1'2  JT  (2i  -  1). 

i=l 

6.4.  Theorem.  If,  for  all  k  e  Z+,  s^  =  t^  +  n^,  where 

tk  df  t(kT)  =  A  sin  (2irfkT  ♦  P)  , 

A,  P,  T  e  R,  0  ^  P  ^  2it,  0  <  f  <  1  /  (2T)  ,  and  nk  df  n(kT)  is  independent,  zero- 
mean,  Gaussian  noise  (nk  has  normal  distribution  with  E(nk)  =  0,  Var(nk)  =  a2, 
and  E(n^n^)  =  E(n^)  E(nj)  =  0  if  i  ?  j),  then  the  following  hold: 

(a)  for  all  k  e  Z+,  sk  has  normal  distribution  with  E(sk)  =  tj.  and 
Var(sk)  =  a2; 

(b)  if  p  =  N  =  2,  then  the  solution  (aj.a  )  of  the  (6.1)  is  given  by  (6.2) 
E(Sj  s^  -  S2  s3)  =  -  A2  sin(2irfT)  sin(4itfT)  , 


and  (6.3)  and 


Var(Sj  s4  -  s2  Sj)  =  2  o4  +  A2  a2  sin2  (2irfkT  +  P)  , 
eCs32  -  s2  s4)  =  a2  +  A2  sin2  (2irfT)  , 

Var (S32  -  s2  s4)  =  3  a4  +  A 2  a2  [sin2  (47rfT  +  P)  +  4  sin2  (6TrfT  ♦  P)  +  sin2  (8irfT  ♦  P)  ] , 
E(s22  -  Sj  S3)  =  a2  +  A2  sin2(2irfT)  , 

and 

Var(s22  -  sx  s3)  =  3  a4  +  a2  a2  [sin2(2TrfT  +  P)  ♦  4  sin2(4irfT  +  P)  +  sin2(6irfT  +  P)] . 

Proof.  The  proof  of  (a)  is  trivial.  We  have  already  shown,  at  the  begin¬ 
ning  of  this  section  that  (6.2)  and  (6.3)  give  the  solution  of  (6.1).  We  now 
complete  the  proof  of  (b). 

Let  NUM1  df  s4  -  s2  s3>  NUM2  df  .  s2  s ^  and  DEN  df  s22  .  ^  Sj.  Now> 
E(NUM1)  =  E((tj  +  nj)(t4  +  n4)  -  (t2  ♦  n2)  (t3  ♦  n3)) 

=  E(tj  t4  -  t2  t3  +  tj  n4  ♦  t4  nx  +  n2  n4  -  t2  n3  -  t3  n£  -  n2  n3) 

*1  t4  ‘  t2  t3  *  tiEfn4^  +  ^Efnj)  +  Efa^)  -  t2E(n3) 
t3  E  (n2)  “  E(n2  n3)  . 

By  hypothesis,  the  expectations  involving  the  ni  are  zero;  thus 

E(NUM1)  =  tj  t4  -  t2  t3. 

Hence,  by  use  of  Corollary  6.2  (with  k  =  m  =  1  and  n  =  2) ,  we  have 

E(NUM1)  =  -  A2  sin(2irfT)  sin(4irfT) . 

We  can  view 

NUM1  -  E(NUM1)  (=  Sj  s4  -  s2  s3  -  (tj  t4  -  t2t3)) 
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as  a  function  of  (s^,  s2>  s^,  s^)  and  expand  in  a  Taylor  series  about 
(tj,  t2,  t3>  t4)  to  obtain 

NUM1  -  E(NUM1)  =  t4  (s1  -  tj)  -  t3  (s2  -  t2)  -  t2(s3  -  tj) 

+  tj  (s4  -  t4)  +  (Sj  -  tj)  (s4  -  t4) 


(s2  "  t2)  (s3  "  *3) 


*4  nl  *3  n2  '  t2  n3  +  *1  n4  +  nl  n4  '  n2  n3’ 


Hence, 


Var(NUMl)  =  E((NUM1  -  E(NUM1))2) 


=  E(t42nI2  *  VV  +  W  +  W  +  W  *  VV 


2t3t4nln2  "  2t2t4nln3  +  2tlt4nln4  +  2  t4  nl2  n4 


2t4nln2n3  +  2t2t3n2n3  “  2tlt3n2n4  '  2  *3  ni  n2  n4 
♦  2t3n22n3  -  2tlt2n3n4  -  2^^^^  ♦  2  t2  n2  n32 


+  2tlnln42  *  2tln2n3n4  -  2nln2n3n4)* 

By  use  of  properties  of  expectation.  Lemma  6.3,  and  hypotheses,  we  have 
Var(NUMl)  =  t^Efaj2)  +  t-32  E(n22)  +  t22E(n32) 


♦  t12E(n42)  +  E(n12n4z)  ♦  E(n22n32) 


=  o2  (tj2  +  t22  +  t32  +  t42)  +  2  a4 


=  2o4  +  A2  a2  j[  sin2(2,nfkT  +  P)  . 
k=l 


1 1 1  - 1 2 


"**  • 1  *  « *  .  *"• »  ■«  •  »T*» 


Next, 


E (NUM2)  =  E((tj  ♦  n3)2  -  (t2  ♦  n2)  (t4  +  n4)) 

E(t32  -  t2t4  +  2t3n3  -  t2n4  -  t4  n2  -  n2  n4  +  n32) 

=  t32  -  t2t4  ♦  2t3E(n3)  -  t2E(n4)  -  t4E(n2)  -  E(n2n4)  +  E(n32) 
=  t32  -  t2t4  +  a2. 

Thus,  by  Corollary  6.2, 

E(NUM2)  =  a2  +  A2  sin2  (2irfT) . 

We  can  view 

NUM2  -  E(NUM2)  (=  s32  -  s2  s4  -  (t32  -  t2  t4  ♦  a2)) 

as  a  function  of  (s2,  s3,  s4)  and  expand  in  a  Taylor  series  about  (t2,  tj,  t4) 
to  obtain 

NUM2  -  E(NUM2)  =  -  a2  -  t4  (*2  -  t2)  +  2  tj  (s3  -  t3)  -  t2  (s4  -  t4) 

(s2  "  ^2)  ^s4  ”  *4)  +  ^s3  ~  ^3) 2 
-a2  -  t4n2  +  2t3n3  -  t2  n4  -  n2  n4  +  n32. 

Proceeding  as  in  the  calculation  of  Var(NUMl),  we  can  show  that 

Var(NUM2)  =  0“  +  t42E(n22)  +  4t32E(n32)  +  t22E(n42)  ♦  E(n22n42) 

+  E(n34)  -  2  a2  E(n32) 

=  a4  +  o2  (t42  +  4t32  +  t22)  +  a4  +  3  a4  -  2  a4 
=  3o4  ♦  a2  (t42  +  4tj2  +  t22) 
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Similarly,  we  can  obtain  expressions  for  E(DEN)  and  Var(DEN)  which  are 
identical  and  almost  identical,  respectively,  to  those  for  E(NUM2)  and 
Var(NUM2).  (Note  that,  by  Corollary  6.2,  t^2  -  t^t^  =  t,2  -  tjt^.) 

As  a  check  on  the  reasonableness  of  the  expectations  in  Theorem  6.4  (b) , 
let  us  consider  some  consequences  of  assuming  that  NUM1,  NUM2 ,  and  DEN 
actually  take  their  expectations  as  values  so  that 

aj(o)  =  (- A2  sin(2nfT)  sin(4nfT))  /  (a2  ♦  A2  sin2  (2irfT) ) 
and 

ai(o)  =  (o‘  +  A2  sin2  (27rfT) )  /  (a2  ♦  A2  sin2  (2nfT) )  “  1. 

(Note:  We  are  not  claiming  that  a^  =  E(NUM1)  /  E(DEN) ,  nor  that  E(aj)  = 

E(NUM1)  /  E(DEN) .)  One  can  readily  show  that 

a.  (0)  =  -2cos(2irfT)  and  a^(0)  =  1, 

which  we  have  shown  |5;  p.  6,  Theorem  6.3  (b)]  gives  the  solution  of  the 
system  (6.1).  Thus,  from  (6.4),  we  get  f  =  f.  Also, 

lim  a^o)  =  0  and  lim  a,(o)  =  1. 

Thus,  in  this  limiting  case,  we  get  f  =  3000  from  (6.4).  (In  this  case,  the 
equation  (2.2)  is  z2  +  1  =  0  which  has  solutions  yielding  0  =  n  /  2  and 
f  =  (6000  /  it)  (it  /  2)  =  3000.)  Therefore,  under  the  above  assumptions,  we  get 
the  same  result,  f  =  f,  as  in  the  (deterministic)  all  signal  -  no  noise  case  if 
o=0  and  we  get  f  =  3000,  which  is  midway  along  the  frequency  band  (0,  6000) 
(and  obviously  minimizes  the  error  f-  f ) ,  in  the  all  noise  -  no  signal  case 
o  . 


The  following  table  gives 
tions  from  Theorem  6.4  (b)  for 


the  expectations,  variances  and  standard  devia- 
f  =  3250,  A  =  1,  and  a  =  0.01,  0.1,  and  1.0. 


a 

1.0 

0.1 

0.01 

E(NUM1) 

0.2566 

0.2566 

0.2566 

Var(NUMl) 

4.1535 

0.02174 

0.0002154 

Std(NUMl) 

2.0380 

0.1473 

0.01468 

E (NUM2) 

1.9830 

0.9930 

0.9831 

Var (NUM2) 

6.7312 

0.03761 

0.0003732 

Std (NUM2) 

2.5945 

0.1939 

0.01932 

E(DEN) 

1.9830 

0.9930 

0.9831 

Var (den) 

3.1045 

0.02114 

0.0002105 

Std (DEN) 

1.7619 

0.1451 

0.01451 

Some  experimenting  with  reasonable  departures  of  NUM1,  NUM2,  and  DEN  from 
their  means  in  case  a  =  1  will  show  that  the  predicted  value  of  f,  f,  given  by 
(6.4),  can  depart  substantially  from  3250. 

Table  1  gives  the  results  of  a  simulation  of  the  signal  and  Gaussian  noise 
case.  The  variables  AM  and  SD  respectively  are  the  arithmetic  mean  and  stan¬ 
dard  deviation  of  100  predicted  values,  f,  of  f.  The  simulation  program  is 
given  as  Appendix  2. 

We  use  a  procedure  given  by  Knuth  [2;  p.  104]  for  generating  random  num¬ 
bers  with  distribution  N(0,  1);  then,  numbers  with  any  normal  distribution  can 
be  generated  by  a  simple  linear  transformation  (if  X  is  N(0,  1),  then  Y  df 
M  +  a  X  is  N(p,o2)).  This  procedure  is  given  below. 


Until  satisfied,  do: 

(1)  Generate  two  uniformly  distributed  random  numbers  r^  and  r^. 

(2)  Set  =  2  Tj  -  1  and  set  V2  =  2  r^-  1. 
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(3) 

(4) 


Set  S  =  Vj2  ♦  V,2. 

If  S  <  1,  set  Xj  ■  (-  2  In  (S)  /  S)^  and  set  x2  = 

L 

V2  (-  2  In  (S)  /  S)  ,  otherwise,  go  to  (1) . 


As  not  all  pairs  of  uniformly  distributed  random  numbers  lead  to  a  pair  of 
normally  distributed  random  numbers,  an  excess  of  such  pairs  must  be  generated. 
The  uniformly  distributed  random  numbers  were  obtained  from  use  of  the  Honey¬ 
well  library  function  RANDT. 


TABLE  1 


Arithmetic  mean  AM  and  standard  deviation  SD  of  100  predicted  values 
f  of  the  Hertz- frequency  f  (=  3250)  of  the  transmitted  signal  in  the 
signal  and  Gaussian  noise  case  for  the  standard  deviation  of  the  noise 
o  =  0.001,  0.01,  0.05,  0.1,  and  0.2  with  the  amplitude  of  the  trans¬ 
mitted  signal  A  =  1,  the  number  of  error  terms  in  the  minimization 
(of  Section  3)  N  =  2,  501,  and  1000,  and  the  number  of  ADC  bits  NB  * 
28.  SNR  is  the  signal-to-noise  ratio  20  log^  (A /a). 


■n 

SNR 

N 

AM 

SD 

0.001 

60.00 

2 

3249.9703 

1.4424 

501 

3249.9993 

0. 5466E-2 

1000 

3249.9995 

0. 2802E-2 

0.01 

40.0 

2 

3249.7444 

14.5075 

501 

3249.9546 

0.5509E-1 

1000 

3249.9389 

0. 2798E-1 

0.05 

26.02 

2 

3249.7308 

72.3807 

501 

3248.9504 

0.2803 

1000 

3248.4615 

0.1399 

0.1 

20.00 

2 

3252.1990 

145.6529 

501 

3245.9330 

0.5945 

1000 

3243.9413 

0.2873 

0.2 

13.98 

2 

3268.8310 

307.4860 

501 

3234.9594 

1.4030 

1000 

3227.3130 

0.6397 
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In  the  no  noise  case  studied  in  [5],  we  showed  that  increasing  N,  the  num¬ 
ber  of  error  terms  in  the  minimization  of  Section  3,  significantly  decreases 
the  variability  of  the  f  values  (see  the  table  on  page  5  of  this  paper).  In 
the  present  simulation,  we  give  N  the  values  2,  501,  and  1000. 

The  program  calls  for  giving  the  standard  deviation  of  the  zero-mean 
Gaussian  noise,  o,  the  sequence  of  values  0.001,  0.01,  0.05,  0.1,  0.2,  0.3, 

0.5.  With  o  *  0.3,  the  variable  DISCR  took  on  some  negative  values  which  led 
to  an  attempt  to  find  the  square  root  of  a  negative  number.  DISCR  has  value 
4.0  A2  -  Al2,  which  is  the  negative  of  the  discriminant  of  the  quadratic 
Z2  +  A1 Z  ♦  A2.  Thus,  with  a  *  0.3,  zeros  of  some  of  the  100  quadratics  are 
real  and  (6.4)  does  not  apply. 

We  note  from  Table  1  that  the  accuracy  of  the  predicted  values,  f,  of  f 
decreases  (3250  -  AM  increases  and  SD  increases)  as  o  increases.  Also,  we  note 
that  the  sequence  of  AM  values  (with  N  =  1000)  3249.9995,  3249.9389,  3248.4615, 
3243.9413,  3227.3130  respectively  corresponding  to  o  =  0.001,  0.01,  0.05,  0.1, 
0.2  are  receding  from  3250  toward  3000  (see  the  comments  on  page  12). 

7.  SIMULATION  IN  THE  SIGNAL  AND  INTERFERENCE  CASE 

Assumo,  for  all  k  e  Z+, 

(7.1)  s(kT)  »  At  sin(2nftkT)  ♦  A^  sin (2nf.kT)  , 

where  f^  is  the  Hertz- frequency  of  the  transmitted  signal  and  f^  is  the  Hertz- 
frequency  of  the  interference.  We  assess,  by  simulation,  the  sampling  varia¬ 
bility  of  the  predicted  values  of  f  and  of  f^  due  to  a  combination  of  quanti¬ 
zation  noise  and  roundoff  noise. 

Suppose  we  wish  ft  to  take  two  distinct  values  (#  F  =  2)  which  arc  50  Hertz 
apart,  say  9975  and  10025,  and  f^  =  10000.  In  accordance  with  the  Sampling 
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Theorem  [ 1 ,  p.  291],  we  must  sample  at  a  rate  of  at  least  20050  (=  2  •  10025) 
samples  per  second.  However,  in  view  of  the  poor  results  (reported  in  [5]) 
obtained  in  case  s(kT)  =  Asin(2nfkT)  and  f  is  "close  to"  an  end-point  of  the 
frequency  band  under  consideration,  we  sampled  at  25000  samples  per  second 
(T  =  1  /  25000) . 

We  give  simulation  results  for  f^  =  10025;  similar  results  were  obtained 
for  f  =  9975.  We  did  not  study  model  behavior  when  f  is  switched  instantane¬ 
ously  from  10025  to  9975  or  vice  versa.  The  simulation  program  is  given  as 
Appendix  3. 

In  case  s(kT)  is  given  by  (7.1),  a  four-pole  model  (2.1)  (p  =  4)  is  appro¬ 
priate  and  the  equation  (2.2)  is  quartic.  Rather  than  use  the  tedious  quartic 
formula  to  solve  (2.2),  we  use  the  Honeywell  library  subprogram  20RP2,  which 
uses  a  modified  Downhill -Newton  method.  For  solving  the  linear  system  (3.2), 
we  use  the  Honeywell  library  subroutine  LINSS,  which  uses  Gauss  elimination 
with  pivoting. 

Both  Z0RP2  and  LINSS  were  modified  as  early  runs  led  to  no  results  ("DEGEN 

ERATE  MATRIX  ..."  error  messages  from  LINSS)  or  poor  results.  Somewhat  better 

results  were  obtained  by  replacing  the  suggested  value,  10  b  or  10  ,  of  the 

-9 

variable  EPS  by  10  .  However,  satisfactory  results  were  obtained  only  upon 

going  to  double-precision  representations  of  data  and  double-precision  arith¬ 
metic  throughout  the  main  program  and  all  the  subprograms  and  decreasing  the 
value  of  EPS  to  10  .  The  Honeywell  version  of  Z0RP2  consists  of  single¬ 

precision  subroutines  DOWNH,  GRAD,  MTALGD,  DIV,  and  POLY.  In  DOWNH,  a  number 

of  logical  IF  statements  involve  constants  10  ^  and  10  ^ .  Upon  going  to  double 

-12  -14 

precision,  these  constants  were  replaced  by  10  and  10  ,  respectively;  how¬ 

ever,  this  resulted  in  "EXP  UFL"  error  messages.  These  constants  were  replaced 
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ultimately  by  10  and  10"  ,  respectively. 

The  program  in  Appendix  3  calculates  f  and  f\  values  from  FC  *  ARCTA (X ,  Y) , 
where  X  and  Y  respectively  are  the  real  part  and  imaginary  part  of  a  zero  of 
the  quartic  and  Y  ^  0.  (For  a  given  complex  conjugate  pair  of  solutions  of 
(2.2),  we  use  the  solution  with  positive  imaginary  part  to  get  a  value  of  f  or 
f^  in  (0,  12500).)  In  case  Y  >_  0,  ARCTA  returns  a  number  in  [0,  it]  so  that 
FC *  ARCTA(X,  Y)  is  in  [0,  12500].  In  case  the  quartic  has  real  zeros  (Y  =  0), 

FC *  ARCTA (X,  Y)  has  value  either  0  or  12500.  The  program  in  Appendix  2  calcu¬ 
lates  f  values  from  FC *  ARCTA(- A1 ,  DSQRT(DISCR) ) ,  where  DISCR  has  as  value  the 
negative  of  the  discriminant  of  the  quadratic  (2.2)  with  p  =  2.  Thus,  as  noted 
in  Section  6,  execution  breaks  down  in  the  case  of  real  zeros  of  the  quadratic 
(for  which  DISCR  has  negative  values). 

We  assume  that  the  amplitudes  (measured  in  volts),  A  and  B  ,  of  the  trans¬ 
mitted  signal  and  the  interference,  respectively,  are  in  [-5,5].  In  the  simu¬ 
lation  of  the  behavior  of  an  NB-bit  analog-to-digital  converter  (=  ADC),  the 
signal  S(J)  in  the  interval  [-5,5],  of  length  10,  is  mapped  into  the  integer 
interval  [-  2^B  *,  2^B  *]  by  following  x  (2^  /  10)  x  with  chopping  to  an 

integer  and  this  integer  is  then  mapped  back  into  a  floating-point  real  number 
in  [-  5,  5]  by  x|-»  (10  /  2NB)  x. 

The  simulation  results  appear  in  Table  2  below.  The  number  of  ADC  bits, 

NB,  takes  values  16,  20,  24,  27,  and  30.  These  values  were  selected  because 
16-bit  ADCs  are  the  most  accurate  ADCs  presently  available  and  20,  24,  27,  and 
30  bits  respectively  correspond  to  approximately  6,  7,  8,  and  9  decimal  digits 
of  precision  (due  to  the  fact  that,  in  a  certain  average  sense,  binary  repre¬ 
sentations  of  floating-point  real  numbers  have  log?(10)  (=  3.32)  times  as  many 
symbols  as  do  the  corresponding  decimal  representations).  The  signal-to- 
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interference  ratio,  SIR,  is  given  by 


SIR  df  lOlogioCA^/Aj2)  =  20  log10(At/A.) 

For  economy.  Table  2  does  not  give  the  simulation  results  for  SIR  <  0,  which 
are  similar  to  those  for  SIR  >  0  except  that  the  approximation  of  f^  (respec¬ 
tively,  f^)  is  better  than  that  of  f^  (respectively,  ft)  for  SIR  >  0  (respec¬ 
tively,  <  0).  More  specifically,  for  large  positive  (respectively,  negative) 
values  of  SIR,  the  arithmetic  mean  of  the  f^  (respectively,  f^)  values  was  much 
closer  to  ft  (respectively,  f^)  than  was  the  arithmetic  mean  of  the  f7  (respec¬ 
tively,  ft)  values  to  fi  (respectively,  ft) .  As  SIR  approaches  0  this  dispar¬ 
ity  diminishes.  This  is  to  be  expected,  as  SIR  >  0  (respectively,  <  0)  implies 
At  >  Ai  (resPectively»  Ai  <  At)-  With  a  27  (or  30)  bit  ADC,  the  predicted 
frequencies,  ft  and  f^ ,  are  quite  accurate  (error  in  the  arithmetic  mean  of  the 
ft  (respectively,  f^)  <  0.00001  (respectively,  0.4))  for  SIR  e  [0,70];  for 
SIR  e  [-70,0],  this  statement  also  holds  if  we  interchange  "f  "  and  "f^". 

This  is  surprising  in  view  of  the  fact  that  for  SIR  =  70  (respectively,  -70), 

At  >  3200  A^  (respectively,  A^  >  3200  At). 

The  program  in  Appendix  3  calls  for  calculating  and  printing  the  arithme¬ 
tic  means  of  the  predicted  magnitudes  of  the  complex  solutions  of  (2.2)  in 
addition  to  the  arithmetic  means  of  the  predicted  frequencies.  The  arithmetic 
means  of  the  magnitudes  of  the  complex  solutions  which  yield  the  predicted 
values  of  f  (=  10025)  are  all  1.00000  to  five  decimal  places  except  for  the 
cases  NB  =  16  and  SIR  =  0,  5,  10,  and  15  respectively  in  which  the  arithmetic 
means  are  0.99978,  0.99974,  0.99985,  and  0.99993.  The  arithmetic  means  of 
the  magnitudes  of  the  complex  solutions  which  yield  the  predicted  values  of 
f^  (=  10000)  are  all  1.00  to  two  decimal  places  except  for  the  cases  NB  =  16 
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and  SIR  =  15,  20,  . . .  ,  80,  NB  «=  20  and  SIR  =  40,  45,  ...  ,80,  and  NB  =  24  and 
SIR  *  65,  70,  75,  80.  These  results  are  expected  as  the  poles  of  the  transfer 
function  should  be  on  the  unit  circle. 


TABLE  2 


Arithmetic  mean  AMp,p  (respectively,  AMpi)  and  standard  deviation 
SDp^,  (respectively,  SDpj)  of  100  predicted  values  of  the  Hertz- 
frequency  of  the  transmitted  signal  (respectively,  interference) 
for  the  signal-to-interference  ratio  SIR  =  0.0,  5.0,  10.0,  ....  80.0 
with  AT  =  5.0,  N  =  496,  and  the  number  of  ADC  bits  NB  =  16,  20,  24, 
27,  and  30. 


NB  =  16 


SIR 

AI 

AMpT 

SDFT 

AMpi 

SDPI 

80.0 

0. 5000E-3 

10025.00001 

0. 14875E-4 

* 

75.0 

0.8891E-3 

10025.00000 

0. 55818E-4 

* 

70.0 

0. 1581E-2 

10025.00001 

0.62049E-4 

* 

65.0 

0. 2812E-2 

10024.99997 

0.62494E-4 

* 

60.0 

0. 5000E-2 

10025.00000 

0. 13758E-3 

* 

55.0 

0.8891E-2 

10024.99997 

0. 21750E-3 

* 

50.0 

0. 1581E-1 

10024.99984 

0.59548E-3 

* 

45.0 

0. 2812E-1 

10024.99931 

0. 11306E-2 

* 

40.0 

0. 5000E-1 

10024.99734 

0. 20244E-2 

* 

35.0 

0. 8891E-1 

10024.99198 

0. 45764E-2 

9904.55117 

0. 12742E+1 

30.0 

0. 1581E+0 

10024.97501 

0. 13919E-1 

9782.96934 

0. 51214E+1 

25.0 

0. 2812E+0 

10024.92350 

0.60974E-1 

9891.21543 

0.20241E+1 

20.0 

0. 5000E+0 

10024.79049 

0. 28133E+0 

9962.73950 

0. 34268E+0 

15.0 

0. 8891E+0 

10024.59744 

0. 93372E+0 

9987.73936 

0. 70680E+0 

10.0 

0. 1581E+1 

10024.53250 

0. 14395E+1 

9996.65213 

0. 12077E+1 

5.0 

0. 2812E+1 

10024.46491 

0.12766E+1 

9998.76546 

0.11762EM 

0.0 

0. 5000E+1 

10024.44931 

0. 70598E+0 

9999.60803 

0. 56777E+0 

AMpi 

is  grossly  in 

error  due  to 

the  occurrence. 

in  calculating 

some  or  all  of 

the  predicted  values  of  FI,  of  real  pairs  (0  and 

12500) ,  rather 

than  complex 

pairs 

,  of  solutions 

of  (2.2). 
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NB  *  20 


SIR 

AI 

SD„_ 

FT 

AMpi 

SDFi 

80.0 

0.5000E-3 

10024.99999 

0. 31531E-5 

* 

75.0 

0.8891E-3 

10025.00000 

0.24865E-4 

* 

70.0 

0. 1581E-2 

10024.99999 

0.60063E-4 

* 

65.0 

0. 2812E-2 

10025.00000 

0. 11830E-3 

* 

60.0 

0. 5000E-2 

10025.00000 

0. 24188E-3 

10299.64857 

0.90224E+1 

55.0 

0.8891E-2 

10024.99994 

0.75772E-3 

9775.32456 

0. 10371E+1 

50.0 

0. 1581E-1 

10024.99960 

0.30767E-2 

9873.41665 

0.10457E+1 

45.0 

0.2812E-1 

10024.99904 

0. 14200E-1 

9955.79372 

0. 18207E+1 

40.0 

0. 5000E-1 

10024.99895 

0.51282E-1 

9985.21355 

0. 87810E+0 

35.0 

0.8891E-1 

10025.00235 

0.88389E-1 

9995.09369 

0.48957E+0 

30.0 

0. 1581E+0 

10025.00380 

0.73353E-1 

9998.37917 

0. 24974E+0 

25.0 

0. 2812E+0 

10024.99760 

0.42521E-1 

9999.50158 

0.74370E-3 

20.0 

0. 5000E+0 

10025.00373 

0. 24424E-1 

9999.85185 

0. 78603E-1 

15.0 

0.8891E+0 

10025.00365 

0. 15603E-1 

9999.94350 

0.40965E-1 

10.0 

0. 1581E+1 

10024.99374 

0.63482E-2 

9999.98351 

0. 88233E-2 

5.0 

0. 2812E+1 

10025.00370 

0.30933E-2 

9999.99406 

0. 13460E- 1 

0.0 

0. 5000E+1 

10025.00169 

0. 11661E-1 

9999.99412 

0.69500E-2 

NB  =  24 

SIR 

A I 

AMpT 

SOpi 

AMpi 

sdfi 

80.0 

0.5000E-3 

10025.00000 

0. 31543E-4 

9813.03648 

0.39144E+1 

75.0 

0. 8891E-3 

10025.00000 

0. 14769E-3 

9853.07626 

0. 14610E+1 

70.0 

0. 1581E-2 

10025.00003 

0.63502E-3 

9943.78711 

0.83820E+0 

65.0 

0. 2812E-2 

10024.99997 

0. 26975E-2 

9981.51822 

0. 83104E+0 

60.0 

0. 5000E-2 

10025.00016 

0. 53971E-2 

9993.55141 

0. 11223E+0 

55.0 

0. 8891E-2 

10025.00026 

0.60723E-2 

9997.80849 

0. 74932E-1 

50.0 

0. 1581E-1 

10025.00035 

0. 31791E-2 

9999.42349 

0. 34453E-1 

45.0 

0. 2812E-1 

10025.00033 

0. 20466E-2 

9999.79271 

0.20292E-1 

40.0 

0. 5000E-1 

10024.99983 

0.19149E-2 

9999.90712 

0. 10978E-2 

35.0 

0.8891E-1 

10025.00032 

0.41177E-3 

9999.97580 

0.67684E-2 

30.0 

0. 1581E+0 

10025.00032 

0. 27273E-3 

9999.99337 

0. 38625E-2 

25.0 

0. 2812E+0 

10025.00032 

0. 17491E-3 

9999.99774 

0. 21642E-2 

20.0 

0. 5000E+0 

10025.00012 

0. 11124E-3 

9999.99949 

0. 46476E-2 

15.0 

0. 8891E+0 

10024.99993 

0. 28941E-3 

10000.00030 

0.36852E-3 

10.0 

0. 1581E+1 

10025.00031 

0.17481E-3 

10000.00010 

0. 39369E-3 

5.0 

0. 2812E+1 

10025.00032 

0.51933E-4 

9999.99995 

0. 221 10E-3 

0.0 

0. 5000E+1 

10025.00001 

0.91795E-3 

9999.99985 

0. 31839E-3 
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NB  =  27 


SIR 

AI 

^FT 

S0FT 

AMpi 

sdfi 

80.0 

0. 5000E-3 

10025.00001 

0. 54492E-3 

9993.09929 

0. 15019E+0 

75.0 

0. 8891E-3 

10025.00004 

0.54903E-3 

9997.66742 

0.38S08E+0 

70.0 

0. 1581E-2 

10025.00003 

0. 24169E-3 

9999.39654 

0. 32871E+0 

65.0 

0.2812E-2 

10025.00000 

0. 28053E-3 

9999.71835 

0. 21289E-1 

60.0 

0. 5000E-2 

10025.00003 

0. 21474E-3 

9999.88056 

0.63278E-1 

55.0 

0. 8891E-2 

10025.00004 

0. 12446E-3 

9999.97034 

0. 19248E-1 

50.0 

0. 1581E-1 

10025.00000 

0. 15211E-3 

9999.97841 

0. 33529E-2 

45.0 

0. 2812E-1 

10024.99999 

0. 41514E-4 

9999.99853 

0. 16168E-1 

40.0 

0. 5000E-1 

10025.00000 

0. 11558E-3 

9999.99448 

0. 10382E-2 

35.0 

0.8891E-1 

10025.00000 

0. 10051E-3 

10000.00254 

0. 57312E-3 

30.0 

0. 1581E+0 

10025.00003 

0. 12463E-3 

9999.99850 

0. 20064E-2 

25.0 

0 . 2812E+0 

10024.99999 

0.28505E-4 

9999.99985 

0. 16137E-2 

20.0 

0. 5000E+0 

10025.00000 

0. 10251E-3 

9999.99950 

0. 10498E-3 

15.0 

0.8891E+0 

10025.00000 

0.95455E-4 

10000.00030 

0.62359E-4 

10.0 

0. 1581E+1 

10025.00000 

0. 11685E-3 

10000.00015 

0.34832E-4 

5.0 

0. 2812E+1 

10025.00002 

0.62149E-4 

9999.99999 

0. 18437E-3 

0.0 

0. 5000E+1 

10024.99998 

0. 24879E-4 

9999.99996 

0. 13870E-4 

NB  =  30 


SIR 

AI 

AMpi 

sdft 

AMFi 

sdfi 

80.0 

0. 5000E-3 

10025.00000 

0. 11038E-2 

10006.70448 

0.49575E-1 

75.0 

0. 8891E-3 

10025.00000 

0. 77445E-3 

10002.07125 

0.86676E-2 

70.0 

0. 1581E-2 

10025.00000 

0.41554E-3 

10000.66534 

0. 13332E-1 

65.0 

0.2812E-2 

10025.00000 

0 . 22524E-3 

10000.20150 

0. 55550E-2 

60.0 

0. 5000E-2 

10025.00000 

0. 12143E-3 

10000.06493 

0. 11716E-2 

55.0 

0. 8891E-2 

10025.00000 

0.68928E-4 

10000.02024 

0. 37293E-3 

50.0 

0. 1581E-1 

10025.00000 

0.49600E-4 

10000.00622 

0. 18744E-3 

45.0 

0. 2812E-1 

10025.00000 

0. 31128E-4 

10000.00295 

0. 22967E-2 

40.0 

0.5000E-1 

10024.99999 

0. 14440E-4 

10000.00057 

0.43003E-3 

35.0 

0. 8891E-1 

10024.99999 

0. 27836E-4 

10000.00057 

0. 32625E-3 

30.0 

0. 1581E+0 

10025.00000 

0. 22309E-4 

10000.00023 

0.40442E-3 

25.0 

0. 2812E+0 

10024.99999 

0. 15418E-4 

10000.00004 

0.80662E-4 

20.0 

0. 5000E+0 

10025.00000 

0. 17013E-4 

10000.00001 

0. 29009E-4 

15.0 

0. 8891E+0 

10025.00000 

0. 10302E-4 

9999.99998 

0.62105E-4 

10.0 

0. 1581E+1 

10024.99999 

0. 13689E-4 

10000.00001 

0.22865E-4 

5.0 

0. 2812E+1 

10025.00000 

0. 25597E-4 

10000.00002 

0. 22428E-4 

0.0 

0. 5000E+1 

10025.00000 

0. 21318E-4 

10000.00001 

0. 19165E-4 
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8.  CONCLUSIONS 


As  indicated  earlier,  the  problem  of  obtaining  confidence  intervals  for 
the  frequencies  of  the  transmitted  signal  appears  to  be  intractable,  both  for 
the  signal  and  noise  case  and  for  the  signal  and  interference  case.  The  case 
of  a  signal  in  noise  and  interference  was  not  considered. 

Simulation  (and  partial  theoretical)  results  in  the  signal  in  independent, 
zero-mean,  Gaussian  noise  case  were  disappointing.  The  writer  is  not  certain 
whether  the  Gaussian  noise  assumption,  so  common  in  the  literature,  is  made 
because  it  is  realistic  or  because  of  its  mathematical  niceties. 

Simulation  in  the  signal  and  interference  case  indicates  the  need  for  an 
ADC  of  accuracy  at  least  27  bits.  At  present,  16-bit  ADCs  are  the  most  accu¬ 
rate  available. 
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APPENDIX  I 


Key  variables  in  the 

Sections  2,  3,  and  6 

TSA  (T) 

FI  (f) 

A  (A) 

NP  (p) 

N  (N) 

NN  (N  ♦  p) 

NB  (NB) 

SIG(K)  (o) 

S(J)  (s  (kT) ) 
P(I,J)  (^(m.N)) 

PO(I)  (*0i(m,N)) 

Al,  A2  (a^,  ap 
DISCR  (4a~-a^2) 

F(K)  (?) 

FC  (  ) 

AM  (AM) 

SD  (SD) 

Key  variables  in  the 

Sections  2,  3,  and  7 

TSA  (T) 

FC  (  ) 

FT  (ft) 

FI  (f4) 


program  in  Appendix  2,  the  corresponding  variables  in 
and  their  interpretations: 


sampling  period  (in  seconds); 

Hertz- frequency  of  the  transmitted  signal; 

amplitude  (in  volts)  of  the  transmitted  signal; 

number  of  poles  of  the  transfer  function; 

number  of  error  terms  in  the  minimization  process 
for  determining  the  prediction  coefficients  aT  in 
(3.1);  3 


number  of  previous  samples  of  the  signal  needed 
to  predict  s (mT) ; 

number  of  ADC  bits; 

standard  deviation  of  the  Gaussian  noise; 


value  of  the  received  signal  at  time  JT; 


coefficient  of 
system  (3.3); 


a. 

3 


in  equation  i  of  the  linear 


constant  on  the  right  in  equation  i  of  the  linear 
system  (3.1); 

coefficients  in  equation  (2.2)  with  p  =  2; 

negative  of  the  discriminant  of  the  quadratic  in 
equation  (2.2)  with  p  =  2; 

predicted  value  of  the  Hertz-frequency  f  of  the 
transmitted  signal; 

radian-to-Hertz  conversion  factor  6000 /it; 
arithmetic  mean  of  100  f  values; 
standard  deviation  of  100  f  values. 


program  in  Appendix  3,  the  corresponding  variables  in 

and  their  interpretations: 

sampling  period  (in  seconds) ; 
radian-to-Hertz  conversion  factor  12500 /w; 

Hertz- frequency  of  the  transmitted  signal; 

Hertz- frequency  of  the  interference; 
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AT 

(At) 

AI 

<Ai> 

NP 

(P) 

N 

(N) 

NN 

(N  +  p) 

NB 

(NB) 

SIR 

(SIR) 

FONE 

FTWO 

(fT) 

amplitude  (in  volts)  of  the  transmitted  signal; 

amplitude  (in  volts)  of  the  interference; 

number  of  poles  of  the  transfer  function; 

number  of  error  terms  in  the  minimization  process 
for  determining  the  prediction  coefficients  ST  in 
(3.1);  J 

number  of  previous  samples  of  the  signal  needed 
to  predict  s(mT); 

number  of  ADC  bits; 

signal-to-interference  ratio  (in  decibels); 

predicted  value  (in  Hertz)  of  the  frequency  of  the 
transmitted  signal; 

predicted  value  (in  Hertz)  of  the  frequency  of  the 
interference. 
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APPENDIX  2 


C  MAXIMUM  ENTROPY  SPECTRAL  DEMODULATION 

C  SIGNAL  AND  (GAUSSIAN)  NOISE  CASE 

C 

DOUBLE  PRECISION  PI , TSA, T PI , A , AM , S D , F 1 , FC , R1  ( I  30 0 ) , R2  ( I  3 00 ) 

DOUBLE  PRECISION  Vl , V 2 , S l , W (  1 300  ) , V l l ( l 300  ) , V 2 l ( I  300 ) , S IG (7 ) 

DOUBLE  PRECISION  X 1  ( 2 6 00  ) , Wl , F J , TC ,TS , U, S ( 1 1  02 ) , PNB , DV, DVI , X 

DOUBLE  PRECISION  P ( 2 , 2 ) , PO ( 2 ) , DEL, DELI , DEL2 , Al , A2 , DISCR ,F ( 1 102) 

COMMON/SAMP/S 

COMMON/PHI/P 

COMMON/ P  HO/ PO 

COMMON/ F  REQ/F 

COMMON/ RAD/ PI 

DATA  PI/ TSA/ 3 . I  4 15926 5 3589 79 32 4D0, 8  3.  3  3 3333  3333  3333 3  3D- 6/ 

DATA  NP/A/AM/SD/Fl/2 , 1 . 0  DO , 0 . 0  DO , 0 . 0D0, 3250.  0  DO/ 

DATA  SIG( i) / SIG ( 2 ) , SIG ( 3 ) , SI G (4 )/0 .0  0ID0, 0 . 0 1D0,0 . 0  5D0,0 . IDO/ 
DATA  SIG (5) , SIG (6) , SIG (7)/0.2D0,0 . 3D0,0. 5D0/ 

TPI=  6. 283185  30  717995  86  48 
FC=  6000. 0D0/PI 

C  GENERATE  GAUSSIAN  NOISE  (0160-0360) 

DO  10  J= 1,1 300 
10  R 1 ( J)=RANDT (1.0) 

DO  20  J= 1,1 300 
20  R2 (J)=RANDT ( 1.0) 

1  =  1 

DO  30  J= 1,1300 

Vl  =  2. 0  D0*R 1  ( J ) -  1 . 0  DO 

V2=2.0D0*R2 (J)-1.0D0 

S1=V1*V 1+V2+V2 

IF (S1.GE.1.0D0)  GO  TO  30 

W (I )=-2. 0  D0*D  LOG (S 1)/S 1 

V11(I)=V1 

V2  l  (I  )  =V2 

1  =  1+  l 

30  CONTINUE 
L-I-l 

LOOP  TO  CALCULATE  AND  PRINT  MEAN  AND  STANDARD 
DEVIATION  OF  100  PREDICTED  FREQUENCIES  FOR 
VARIOUS  VALUES  OF  SIGMA 
DO  4  0  K  1=1,7 
DO  50  J=l,L 


o  o  o  no 


Z=SIG(K l)*DSQRT  (W(  J)  ) 

XL  (  2*  J  -  1 )  =V  1 1 ( J) *  Z 
50  X  1  ( 2*  J ) =V  2 1  ( J ) *  Z 

C  CALCULATE  TRANSMITTED  SIGNAL  PLUS  NOISE 

Wl=F  l*T  PI 
DO  6  0  J«  1,1102 
FJ=J 

TC=FJ*TSA 
TS=W  I  *TC 
U=A*DSIN(TS) 

60  S(J)=U+X1(J) 

C  SIMULATE  ANALOG-TO-DIGITAL  CONVERSION  OF  SIGNAL 

DO  45  NB=  16,28,  12 
WRITE<6,1)  NB 

1  FORMAT(  IX,  "NUMBER  OF  BITS  =  *,l5/) 

N  BP=  2*  *NB 

PNB=NBP 
DV=  10. 0D0/PNB 
DVI  =  PNB/ I  0 . 0  DO 
DO  70  J=l,  1102 
X=S  < J) 

X=X*DVI 
K  P=  X 
X=KP 

70  S(J)=X*DV 

LOOP  TO  VARY  N  <=  THE  NUMBER  OF  ERROR  TERMS  IN 
THE  MINIMIZATION  OF  SECTION  3) 

DO  55  N=2, 1000,499 
N  N=N  P+N 
KL=NN+  I 
KU=KL+  99 

LOOP  TO  GENERATE  100  (=  KU-KL+1)  PREDICTED  FREQUENCIES 
DO  80  K=KL, KU 
KK  =  K 

CALCULATE  COEFFICIENTS  FOR  LINEAR  SYSTEM  (3.2) 

USING  SUBROUTINE  COV  (COVARIANCE  METHOD) 

CALL  COV  <NP,NN,KK) 

C  SOLVE  LINEAR  SYSTEM  (3.2)  BY  CRAMER'S  RULE 

DEL=P(1, 1)*P  (2,2)-P(l,2)*P  (2,1) 

DELl  =  P(l,2  )*PO(2)~P  (2,2  )*PO(l) 

DEL2  =  P  (2,  1  )*PO  (1) -P  (  1,  l )  * PO  ( 2 ) 

A l  =  DEL  l/D  EL 
A2=DEL2/DEL 

C  SOLVE  EQUATION  (2.2)  BY  QUADRATIC  FORMULA 

DISCR=  4.0D0*A2-A1*A1 
I F (DISCR.GE  .0 . 0  DO )  GO  TO  80 
WR ITE (6,2)  K 

2  FORMAT ( IX, "DISCR  NEG  AT  STEP  ",I8) 

CALCULATE  PREDICTED  FREQUENCY  BY  USE  OF  (6.4) 

•OF (K )-FC*ARCTA(“Al, DSQRT (DISCR) ) 
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CALCULATE  AND  PRINT  MEAN  AND  STANDARD  DEVIATION 
OF  PREDICTED  FREQUENCIES 
CALL  STATS  (AM,  SD,KL,KU) 

WRITE<6,3)  N, AM,SP 
3  FORMAT ( IX, l5,2G24.8/) 

55  CONTINUE 
45  CONTINUE 
40  CONTINUE 
STOP 
END 
C 

SUBROUTINE  COV ( NP , NN , LP) 

C 

DOUBLE  PRECISION  S  (  I  102 )  , P  (  2 , 2 ) , PO (2 ) , B 

COMMON/SAMP/S 

COMMON/ PHI/P 

COMMON/ P  HQ/ PO 

L-LP-  1 

NI-NN-NP 

NL-LP-NI 

B-0. 0D0 

C  LOOPS  TO  CALCULATE  PHI(J,J),  1<»J<-NP 

DO  10  J«NL,L 

10  B-B+S ( J)*S ( J) 

DO  11  J-  l , NP 

K-LP-J 

1-NL-J 

B-B+S  (I)*S  < I ) -S  (K)*S  (K) 

11  P  (J,  J  )-B 

DO  12  KK-l,NP 
B- 0.0  DO 

C  LOOP  TO  CALCULATE  PHI(0,KK),  1<-KK<»NP  /  STORE  IN  PO(KK) 

DO  13  J-1,NI 
N-LP-J 
M-N-KK 

13  B-B+S (N)*S (M) 

PO  (KK)  =  B 

C  LOOP  TO  CALCULATE  PHI(I,K)  FOR  I  NOT  =  K 

I F (KK. EQ. NP)  GO  TO  12 
DO  14  J-  l , NP-KK 
I  “J 

K-KK+J 
N-LP-J 
M-N-KK 
Nl-NL-J 
Ml -N  l-KK 

B-B+S (Nl)*S (Ml)-S (N)*S (M) 

P  (I  /  K  )  -B 

14  P  (K, I ) - B 
12  CONTINUE 

RETURN 

END 
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SUBROUTINE  STATS  (AM,SD,KL,KU) 

STATS  CALCULATES  THE  MEAN  AND  STANDARD  DEVIATION 
OF  KU-KL+i  PREDICTED  FREQUENCIES  F(J) 

DOUBLE  PRECISION  F ( l 10 2 ) , S l , S 2 , RM , AM ,SD 
COMMON/F  REQ/F 
S 1=0. 0D0 
S2  =  0. 0D0 
RM=KU-KL+I 
DO  10  J=KL,KU 
10  S  1=S  l+F  ( J) 

AM=S 1/RM 
DO  20  J=KL,KU 

20  S2=S2+  (F  (J)-AM)*  (F  (J)-AM) 

SD=DSQRT  (S2/RM) 

RETURN 

END 

FUNCTION  ARCTA ( X , Y ) 

SUBPROGRAM  TO  CALCULATE  ARCTANGENT 
VALUES  IN  -PI  TO  PI 
DOUBLE  PRECISION  PI , HPI , X , Y , ARCTA 
COMMON/RAD/PI 
HPI= 1.57079632 67 9489662D0 
I F  ( X )  1,2,3 

1  A  RCTA=  DAT  AN  (  Y/X  )  +  PI 
RETURN 

2  ARCTA* HPI 
RETURN 

3  ARCTA-DATAN (Y/X) 

RETURN 

END 


; 


1 1 1 - 30 


.  ..  MU 


APPENDIX  3 


C  MAXIMUM  ENTROPY  SPECTRAL  DEMODULATION 

C  SIGNAL  AND  INTERFERENCE  CASE 

C 

DOUBLE  PRECISION  PI , TPI ,TSA , FC , AM, SD , AT , AI , FT , F I , S IR , WT ,WI 
DOUBLE  PRECISION  F J , TC , TST , TS I , S ( 5 0 0 ) , PNB , D V , DVI , X , Y , A ( l 0 ) 
DOUBLE  PRECISION  P ( 4 , 4 ) , PO ( 4 ) , RR ( 5 ) , C R ( 5 ) , B ( 5 ) , A UX  (  4 ) , R ( 4 ) 
DOUBLE  PRECISION  T EMP , RMAG 1 ( 1 0 0 ) , RMAG2 ( 1 0 0 ) , FONE ( 1 0 0 ) 
DOUBLE  PRECISION  FT WO (  10 0 ) , F 1 l ( 2 ) , Rl 1 ( 2 ) 

COMMON/SAMP/S 
COMMON/PHI/P 
COMMON/ P HO/ PO 
COMMON/F  REQ/FONE, FTWO 
COMMON/MAG/RMAG  l,  RMAG2 
COMMON/S IME/A ,R,AUX 
COMMON/ BRC/B,RR,CR 

DATA  AM, SD, AT, FT, FI/ 0.0  DO, 0.0  DO ,5.0  DO, 10025.0  DO,  10  000. 0  DO/ 
DATA  NP  ,  N/  4 , 4  9  6/ 

PI  =  3.  1  4  15926  5358979  324 
TP  1=6. 2  8  3185  30  7  L7  958  6  47 
TSA=  l .  0  DO/ 2  5  000 . 0  DO 
FC=  12  50  0. 0  D0/PI 
NN=N  P+N 
KL=NN+ l 

C  LOOP  TO  GIVE  SIR  VALUES  60,55,50,  ...  ,5,0 

S  IR=  6  0. 0D0 
1 l  WRITE (6,201)  SIR 
201  FORMAT( IX, "SIR  =",Gl7.8//) 

AI=AT*10.0D0**(-SIR/20.0D0) 

W  RITE ( 6 , 1 11)  AI 
111  FORMAT(  IX,  "AI  =",G17.8//) 

WT  =  FT*T  PI 
W I  =  F I *T  PI 

C  LOOP  TO  CALCULATE  100  PREDICTED  VALUES  OF  FT  AND  OF  FI 

DO  33  K=  l ,  100 
DO  3  4  J=  1,5  00 
F  J  =  J  +  50  0*  (K-  1) 

TC=FJ*TSA 

TST=WT*TC 

TSI=WI*TC 

C  CALCULATE  RECEIVED  SIGNAL  FROM  (7.1) 

S ( J ) =AT*  DS I N (TST ) +A I *  DSI N (TSI) 

34  CONTINUE 
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C  SIMULATE  16-BIT  ANALOG-TO-DIGITAL  CONVERSION  OF  SIGNAL 

N  B=  16 
NBP=2*  *NB 
PN  B=N  BP 
DV=  LO.  ODO/PNB 
DVI=PNB/ 10. 0D0 
DO  15  J=  1,500 
X=S (J) 

X=X*DVI 
K  P=  X 
X  =  KP 

15  S ( J)=X*DV 

CALCULATE  COEFFICIENTS  FOR  LINEAR  SYSTEM  (3.2) 

USING  SUBROUTINE  COV  (COVARIANCE  METHOD) 

CALL  COV (NP,NN,KL) 

L=0 

DO  L  J=l,NP 
DO  9  1=  L,  J 
LL=L+I 

A (LL) =P (I , J) 

9  CONTINUE 
L=L+J 

R  ( J)  =  -PO  (J) 

1  CONTINUE 
I  E  R=  0 

C  SOLVE  LINEAR  SYSTEM  (3.2)  USING  HONEYWELL  SUBROUTINE  LINSS 

CALL  LINSS (NP, IER) 

IF(IER)  2,3,5 
5  NQ=NP-IER 

WR ITE  (  6,202)  NQ , IER 

202  FORMAT ( IX,  "DEGENERATE  MATRIX  OF  DE GE NE RACY " , I  8 ,  "  RANK  =",I8) 
GO  TO  33 

2  WRITE  (  6,2  0  3) 

20  3  FORMAT  (  IX,  "MATRI  X  POSSIBLY  SINGULAR") 

GO  TO  33 

C  ASSIGN  ELEMENTS  IN  SOLUTION  4-TUPLE  OF  SYSTEM  (3.2) 

C  (DENOTED  BY  "R(J)M  IN  LINSS  IN  PLACE  OF  "A(J)"  IN  (3.2)) 

C  TO  B ( J )  WITH  B( 5)=1, B( 4 )=R ( 1) , B( 3) =R (2) , B( 2 )=R ( 3) , 

C  B  (  l  )  =R  (  4  )  AND  SOLVE  THE  EQUATION  B(5)*Z**4  +  B(4)*Z**3 

C  +  B (  3  )  *  Z  *  *2  +  B ( 2 ) *  Z  +  B(l)  =  0  OF  THE  FORM  (2.2) 

C  BY  USE  OF  HONEYWELL  SUBPROGRAM  ZORP2  WHICH  CONSISTS 

C  OF  SUBROUTINES  DOWNH , G RAD , MT ALG D, DI V ,  AND  POLY 

3  NAC=NP+l 

B ( NAC)  =  1 . 0  DO 
DO  7  J=  L , NP 
7  B ( J)=R (NAC-J) 

CALL  DOWNH ( B ,NP, RR,CR) 

1=1 

DO  44  L=  1, NP 

IF  (CR(  L)  .  LT  .0 . 0D0)  GO  TO  44 
X  =  R  R ( L) 

Y=C  R  (  L) 

R l  1  (I ) =DSQRT (X*X  +  Y*Y) 
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C  CALCULATE  PREDICTED  FREQUENCY 

Fl 1(1) =FC*A  RCTA ( X, Y ) 

1  =  1+  L 

4  4  CONTINUE 

IF(FLl  (1)  .LE.Fl  L  (2)  )  GO  TO  4  3 
TEMP=F  11(1) 

Fill  l)=Fll(2) 

Fl  1  (2)  =TEMP 
TEMP=R  11(1) 

R  1  1  (  1 )  =R  11(2) 

Rl  l  (2  )  =T EMP 
43  RMAG1  (  K  )  =R  1  1  (  l) 

RMAG2 ( K ) =R 1 1 (2) 

FONE  (  K  )  =F  1  1  (1) 

FT  WO (K ) =F 1  1  (  2) 

33  CONTINUE 

CALCULATE  AND  PRINT  MEAN  AND  STANDARD  DEVIATION 
OF  100  PREDICTED  VALUES  OF  RMAG  l , RMAG2 , F 1 ,  AND  F2 
CALL  STATS (AM, SD, RMAGl) 

WRITE  (  6,4  00)  AM,SD 
400  FORMAT  (  IX,  2D25.  18) 

CALL  STATS  (AM, SD,RMAG2) 

WR IT  E (  6 , 4  0  0  )  AM , S  D 
CALL  STATS (AM, SD, FONE) 

WRITE  (6,4  00)  AM  ,SD 
CALL  STATS  (AM, SD,FTWO) 

WR ITE(6,400)  AM , SD 
WRITE  (6,8  88) 

8  88  FORMAT  (  IX,//) 

SIR=SIR-5.0D0 
IF (SIR.LT.0.0D0)  STOP 
GO  TO  11 
END 
C 

SUBROUTINE  COV  (  NP  ,  NN  ,  LP) 

C 

C  COV  WAS  WRITTEN  BY  CAPT .  KENNETH  WILSON  OF  RADC  TO 

C  IMPLEMENT  THE  COVARIANCE  METHOD  OF  LINEAR  PREDICTION 

C  (SEE  MAKHOUL  (2,  P.  564)) 

DOUBLE  PRECISION  S ( 50 0 ) , P ( 4 , 4 ) , PO ( 4 ) , A 

COMMON/SAMP/S 

COMMON/PHI/P 

COMMON/PHO/PO 

L=LP-1 

NI=NN-NP 

NL=LP-NI 

A=  0 . 0  DO 

C  LOOPS  TO  CALCULATE  PHI(J,J),  1<=J<=NP 

DO  10  J=NL,L 

10  A=A+S ( J)*S ( J) 

DO  11  J=  1,  NP 

K=LP-J 

I=NL-J 

A=A  +  S  (I)+S  (I)-S ( K ) * S  (K) 

11  P(J,J)=A 
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DO  12  KK=1,NP 
A=  0 . 0  DO 

LOOP  TO  CALCULATE  PHI(0,KK),  1<KK<NP  /  STORE  IN  PO(KK) 
DO  13  J=l,NI 
N=LP-J 
M=N-KK 

1 3  A~A+S (N)*S (M) 

PO (KK) =A 

LOOP  TO  CALCULATE  PHI(I,K)  FOR  I  NOT  =  K 
IF(KK.EQ.NP)  GO  TO  12 
DO  14  J=1,NP-KK 
I=J 

K=KK+J 

N=LP-J 

M=N-KK 

N1=NL-J 

M1=N1"KK 

A=A+S (Nl)*S(Ml)-S(N)*S(M) 

P (I  ,K)=A 

14  P ( K, I  )  =A 
12  CONTINUE 

RETURN 

END 


SUBROUTINE  LINSS(M,IER) 


LINSS  IS  A  DOUBLE  PRECISION  VERSION  OF  THE  HONEYWELL 
SUBROUTINE  LINSS  FOR  SOLVING  A  LINEAR  SYSTEM  WITH 
SYMMETRIC  COEFFICIENT  MATRIX.  LINSS  USES  GAUSS  ELIMINATION 
WITH  PIVOTING  IN  THE  MAIN  DIAGONAL  ONLY,  TO  PRESERVE 
SYMMETRY.  IER  IS  AN  ERROR  RETURN  AS  FOLLOWS:  IER=0 
INDICATES  NO  ERROR;  IER=-1  INDICATES  NO  RESULT  AS  NP< 0 
OR  A  PIVOT  ELEMENT  WAS  ZERO  DURING  ELIMINATION;  IER=K 
IS  A  WARNING  OF  POSSIBLE  LOSS  OF  SIGNIFICANCE  (OF  L 
SIGNIFICANT  DIGITS  IF  EPS  -  10**(“L))  AT  ELIMINATION 

STEP  K+l  AND  ,  WITH  WELL-CONDITIONED  A  AND  APPROPRIATE 
EPS,  THAT  A  MAY  HAVE  A  RANK  OF  K. 

DOUBLE  PRECISION  A ( 1 0 ) , R ( 4 ) , AUX ( 4 ) , PI V , TB , TOL , PI VI , EPS 
COMMON/S  IME/A  ,  R,  AUX 
E  PS=  1.0D-  18 
IF(M.LE.O)  GO  TO  24 
SEARCH  FOR  PIVOT 
I  ER=  0 
PIV=0. ODO 
L=0 

DO  3  K= 1 , M 
L=L+K 

T  B=  DABS (A (L) ) 

IF(TB-PIV)  3,3,2 
PIV=TB 
I  =  L 
J=K 

CONTINUE 
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TOL=E  PS*  P  I V 
LST=  0 
LEND=M-1 

C  ELIMINATION  LOOP 

DO  18  K=  1 ,  M 
IF(PIV)  24,24,4 

4  IF(IER)  7,5,7 

5  IF(PIV-TOL)  6,6,7 

6  IER=K-I 

7  LT= J-K 
LST=  LST+  K 

C  PIVOT  ROW  REDUCTION  AND  ROW  INTERCHANGE  IN  RIGHT  SIDE  R 

PIVI= I.ODO/A (I) 

TB=PIVI*R  ( J) 

R(J)=R(K) 

R(K)=TB 

IF(K-M)  9,19,19 

ROW  AND  COLUMN  INTERCHANGE  AND  PIVOT  ROW  REDUCTION  A 
PIVOT  COLUMN  SAVED  IN  AUX 
9  LR=LST+  (LT*  (K+J-l)  )/2 
LL=LR 
L=LST 

DO  14  I  1=  K , LEND 
L=L+ I I 
LL=LL+1 

IF(L-LR)  12,10,11 

10  A (LL)=A ( LST) 

TB=A  (L) 

GO  TO  13 

11  LL=L+LT 

12  TB=A  (  LL) 

A (LL) =A (L) 

13  AUX  (II )  =TB 

14  A ( L) =PI VI*TB 
A ( LST ) =LT 

C  ELEMENT  REDUCTION  AND  SEARCH  FOR  NEXT  PIVOT 

PI  V=0. 0D0 
LLST=LST 
LT=  0 

DO  18  I  1=  K , LEND 
P I VI  = -AUX  (II) 

LL=LLST 
LT=  LT+  l 

DO  15  LLD=I I , LEND 

LL=LL+LLD 

L=LL+LT 

15  A (L)=A (L)+PIVI*A (LL) 

LLST=LLST+I I 
LR=LLST+LT 

T  B=  DABS ( A ( LR)  ) 

IF(TB-PIV)  17,17,16 

16  PI  V=TB 
I  =  LR 


J=I  1+1 

17  LL= 1 1+ 1 

18  R(LL)=R(LL)+PIVI*R(K) 

C  BACK  SOLUTION  AND  INTERCHANGE 

19  I F ( LEND )  24,23,20 
2  0  I  I=M 

DO  22  1  =  2, M 
LST=LST-I  I 
I  1=11-1 

L=A (LST) +0. 5  DO 
TB=R (II) 

LL=  1 1 
K=LST 

DO  21  LT= II, LEND 
L  L=  L  L+  1 


K=K+LT 

TB=TB-A  (  K  )  *  R  (LL) 

K=I  I+L 
R  ( 1 1 )  =R ( K) 

R  (K)=TB 
RETURN 
I E  R=  - 1 
RETURN 
END 

SUBROUTINE  DOWNH ( A , NAR, RR , CR ) 

DOUBLE  PRECISION  A  (  10  )  ,  RR  (  5)  ,  CR  (  5  )  ,  Q  (  10  )  ,  B  (  3) 

DOUBLE  PRECISION  ANP P, DI SC , X , Y , C 

J  =  0 

N  =  NAR 

NPLl=N+l 

AN  PP=A (NPL1) 

DO  102  I =  l , N PLl 

IF  (A  ( I )  )  10  3,  102,  103 

CONTINUE 

C=  DABS (A ( I )/A ( NPLl ) ) 

LU=  120 
LL=-  120 

IF(C-2.0D0**LU)  100,  100,  101 

IF(C-2.0D0**LL)  101,105,105 

NAR=  -NAR 
GO  TO  5001 
I  1=  (LU+LL)/2 

IF (C-2.0D0**II)  110,110,109 

L  L=  1 1 
GO  TO  1 11 
LU=I  I 

IF(LU-LL-l)  500  1,  1  12,  105 
I  B=  1 1  /  N 

IF (IB)  114,120,114 
DO  115  1=1, NPLl 
11=1-1 

A(I)  =A  (I)*(2.0  D0**(H*IB)) 
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120  DO  121  J  1=  l ,  NPLl 
12  1  A  (  Jl)  =A  ( J1  )/A  (NPLl) 

201  I F  (N )  2001,  2001,206 

206  I  F  ( A  (  1 )  )  30  l,  2  11,  301 

211  J=J+1 

RR(J) =0. 0D0 
CR(  J)  =0. 0D0 
DO  22  1  J  L=  1,  N 
22  1  A  ( J1)=A  (Jl+1) 

N  =  N-  1 
GO  TO  201 

30  L  IF  (  N~  2 )  60  1,5  0  1,4  01 

401  CALL  GRAD (A,N,X, Y) 

421  IF  (DABS  (Y)-DABS  (X*  1.  OD-8  )  )  4  31,4  3  1,4  4  1 

431  Y=  0.0  DO 
44  1  J  =  J+1 

RR(J)=X 
C  R ( J) =Y 

I  F  (  Y  )  4  6  1,  102  1,461 
461  J=J+1 

RR( J)=X 
CR( J)=-Y 
GO  TO  10  11 

501  DISC=A (2 )* *2-4. 0D0*A (  1) 

IF ( DISC)  521,541,541 
521  Y=DSQRT(-DISC)/2.0D0 
X=-A  (  2  )/2 . 0D0 
GO  TO  421 
54  1  J=J+1 

RR(J)  =  (-A(2)+DSQRT(DISC)  )  /  2 . 0  DO 
C  R  (  J  )  =  0 .  0D0 
GO  TO  1021 
60  1  J=  J+l 

RR(J)=-A  (  1) 

CK(  J)=0.  ODO 
GO  TO  2001 
10  ll  B(  l)  =  X**2+Y**2 
B ( 2)=-2. 0D0*X 
B  (  3  )  =  l . 0  DO 
N  B=  2 

GO  TO  1041 
1021  B (  1 ) =-RR ( J) 

B  (  2  )  =  1 . 0  DO 
NB=  1 

1041  CALL  DIV (A, B, N,NB,Q) 

DO  106  1  Jl= 1,N 
106  I  A (J1)=Q( Jl) 

IF  (C  R  (  J )  )  108  l,  107  1,  1081 

1071  N=N-l 

GO  TO  201 
10  8  l  N=N- 2 

GO  TO  201 

2001  IF  (IB)  200  2,2005,2002 
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200  2  DO  2  0  00  1=1,  NAR 

RR(I)=RR(D*  (2.0  DO**  (IB)  ) 

2000  CR(I)=CR(I)* (2.0D0**(IB)  ) 

2005  NPl=NAR+l 

DO  2011  1=2, NP1 
2011  A  ( I )  =  0 . 0  DO 
A  (  l)  =  1 . 0  DO 
NA=  0 
J  =  1 

2021  IF(CR(J>)  204  1,2061,2041 
2041  NB=  2 

B  (  3)  =  1.0  DO 
B(2)=-2.0D0*RR(J) 

B(  l)  =  RR(J)**2+CR(J)**2 
J  =  J+2 

GO  TO  2081 
2061  NB=1 

B ( 2)= 1. 0D0 
B  (  L)=-RR(  J)  . 

J=J+1 

2081  CALL  MTALGD (A , NA , B , NB ,Q ) 

NA=N  B+NA 
NAPLl=NA+ 1 
DO  2091  1=  1,  NA PLl 
2091  A  ( I ) -Q ( I ) 

IF(NA-NAR)  202  1,  3001,  3001 
300  1  DO  3011  J2=  l , NPLl 
3011  A( J2)=A (J2)*ANPP 
5001  RETURN 
END 

SUBROUTINE  GRAD  (A  ,  N  ,  X  Z  ,  Y  Z  ) 

DOUBLE  PRECISION  A (  1  0 )  , X ( 3 ) , Y ( 3 ) , RP ( 3 ) , C P ( 3 ) , RHO ( 3 ) , PHI  ( 3 ) 
DOUBLE  PRECISION  ABS P ( 3 ) , PR ( 3 ) , PC { 3 ) , P I , X 2 , Y Z , RHO Z , PHI Z , S U , U 
DOUBLE  PRECISION  PS  I , TO P , BOT , COS  I , S IN E , DZ , ABS P Z , PR Z , PC Z , RZ 
DOUBLE  PRECISION  C Z , THETA , DTHETA , R HOS , PHIS 
PI=  3.  14  1592  6  535  8  979  324 
MTST=  1 
10  l  XZ=  0. 0D0 
YZ=  1 . 0D0 
DZ=  2 . 0  DO 
RHOZ=  1. 0D0 
PHI  Z=PI/2 . 0D0 

201  CALL  POLY (N,A,RZ ,CZ, PRZ,PCZ , RHOZ , PH IZ ) 

221  S  U=DSQRT ( PRZ*  *2  +  PCZ*  *2 ) 

A  BSP  Z=  DSQRT  (RZ**2+CZ**2) 

U=2.0D0*ABSPZ*SU 

PSI=DATAN(U) 

TOP=RZ*PCZ-CZ*PRZ 
BOT=-  (R  Z*  P  RZ+C  Z*PCZ) 

TH  ETA=A RCTA ( BOT , TOP) 

COS 1=  DCOS (THETA+PHIZ) 
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SINE=DSIN  (THETA+PHIZ) 
JF(ABSPZ)  300,5001,300 


300  IF(SU)  301,50  1,  301 
30  1  IF(RHOZ)  32  1,4  0  1,  321 

321  IF  (ABSPZ/  ( RHOZ*S  U)  “  l.  0  D-  8  )  500  1,5001,  7  01 
351  IF  (ABSPZ/ (  RHOZ*S  U)  -  10. 0D0*  *  (-MTST)  )  80  1,8  0  1,4  01 
401  DZ=  DZ/  8 . 0  DO 
I  M=  0 

DO  431  1=1,3 
DZ=2.0D0*DZ 
X(I  )=XZ+DZ*COSI 


Y  (I  )=YZ+DZ*SINE 

RHO (I ) =DSQRT (X(I)**2+Y(I)**2) 

P  H I  (  I  )  =A  RCTA  (  X  ( I  )  ,  Y  ( I  )  ) 

CALL  POLY  (N  ,A,RP(  I  )  ,CP(I  )  ,  PR  (I  )  ,  PC  (I )  ,  RHO  (I  )  ,  PHI  (I  )  ) 
A  BS  P ( I ) =DSQRT ( RP ( I ) *  *2  +C  P ( I ) *  *2 ) 

IF  (ABS  PZ-ABSP  ( I  )  )  4  31,4  31,4  21 

42  1  ABSPZ  =  ABSP(I) 

I  M=  I 

43  1  CONT INUE 

IF(IM)  44  1,4  4  1,4  61 
441  DZ=  DZ/ 8 . 0  DO 

IF(RHOZ)  44  3,4  45,4  43 

44  3  IF(DZ/RHOZ- l.OD-8)  4  51,4  5  1,4  01 

445  IF  (DZ-  1.0D-8)  45  1,4  5  1,4  0  1 

451  IF(SU-ABSPZ)  50  1,50  1,5001 
461  DZ=  (2.0D0**(  IM-2)  )*DZ 
XZ=X  (I M) 

YZ= Y (IM) 

PHIZ=PHI (IM) 

PRZ=PR  (  IM) 

PCZ=PC  (  IM) 

RHOZ=RHO  (IM) 

R  Z  =  R  P ( I M) 

C  Z  =C  P  ( I  M) 

GOTO  221 
501  DZ=  1 . 0  DO 

DT  HETA=  PI/  10. 0D0 
521  T  HETA=  0 . 0  DO 

DO  561  1=  1,20 
T  HETA=T  HETA+  DTHETA 
XS=XZ+DZ*DCOS  (PHIZ+THETA) 

YS=YZ+DZ*DSIN(PHI  Z+THETA) 

RHOS=DSQRT  (XS*  *2+YS*  *2) 

PHI S= A  RCT A ( XS , YS ) 

CALL  POLY  (N,A  ,  RS  ,CS  ,  PRS,  PCS  ,  RHOS  ,  PHIS) 

A  BS  P ( l )=DSQRT (RS**2+CS **2) 

IF  (ABSPZ-ABSP(  1  )  )  56  1,56  1,6  0  1 
56  1  CONTINUE 

DZ=  DZ/2 .0D0 
IF (RHOS)  563,565,563 
56  3  I  F  (DZ/RHOS- 1.  OD-8)  5  0  0  1,500  1,521 


YZ=  YS 
PHIZ- PHIS 
RHOZ=R  HOS 
ABSPZ=ABSP(l  ) 

PRZ=PRS 
PCZ=  PCS 
RZ=RS 
CZ=CS 
GO  TO  221 

701  IF  (PSI- l.OD-7)  711,711,351 
7  L  1  IF(SU-ABSPZ)  50  1,50  1,  351 
80  1  RHO  (  1)  =RHOZ+BOT/S  U*  *2 
IF  (RHO(  1 )  )  90  1,9  0  1,  8  16 
816  PHI(l)=PHI 2+TOP/ ( RHOZ*  S  U*  *  2 ) 

82  1  CALL  POLY (N ,A,RZ ,CZ ,PRZ,PCZ ,RHO ( 1)  , PHI  ( 1) ) 

ABS  P (  1)=DSQRT(RZ**2+CZ**2) 

IF  (ABSP  (  1  ) -ABSPZ)  8  51,8  81,881 
841  XZ=RHOZ*DCOS (PHIZ) 

YZ=RHOZ*DSIN(  PHIZ) 

GO  TO  5001 
851  RHOZ  =  RHO(l) 

A  BS  PZ  =  A  BS  P  (  1) 

PHI  Z=  PHI  (1) 

TOP=RZ*PCZ-CZ*PRZ 
BOT=-  (R  Z*  PRZ+C  Z*  PCZ ) 

S  U=  DSQRT ( PRZ*  *2  + PCZ*  *2) 

IF(SU)  855,50  1,8  55 
855  U=2 . 0D0*ABSPZ*SU 
PS I  =  DAT AN ( U ) 

IF  (ABSPZ/  (RHOZ*SU) -10. 0D0**  (-MTST)  )  861,861,901 
86  1  IF  (ABSPZ/ (RHOZ*S  U) -1.  OD-8  )  8  4  1,8  4  1,8  /1 
871  IF(PSI- l.OD-7)  88  1,8  8  1,  8  01 

88  1  IF  (SU- ABSPZ)  501,501,901 
90  1  D  Z=A  BS  PZ/S  U 

XZ=RHOZ*DCOS (PHIZ) 

YZ=RHOZ*DSIN (PHIZ) 

MTST=MTST+  l 
GO  TO  201 
5  00  1  RETURN 
END 

S  UBROUTINE  MTA LGD (AARG , NA, BARG, NB,C) 

DOUBLE  PRECISION  AARG ( 1 0 ), BARG ( 1 0 ), C ( 1 0 ), A ( 1 0 ), B ( 1 0 ), TEMP 
1  NAPL1-NA+ 1 

DO  21  Jl= 1, NAPLl 
21  A  (  J  1)=AARG  (  Jl> 

NBPH=N  B+  1 
DO  41  Jl= l, NBPLl 
41  B(Jl)  =  BARG(Jl) 

NCPL1=NAPL1+NBPL1-1 
DO  91  Jl»l,NCPLl 


TEMP=0.  ODO 
DO  81  J2  = l , J 1 
IF1J2-NAPL1)  61,61,81 
N  2  =  J  1- J2+  1 

IF(N2-NBPH)  71,71,81 
TEMP=T  EMP+A (J2)*B(N2) 

CONTINUE 
C ( J1 ) =TEMP 
CONTINUE 
RETURN 
END 

SUBROUTINE  DI V  ( A , B , N A , NB , Q ) 

DOUBLE  PRECISION  A ( 1 0 ) , B ( l 0 ) , Q ( 1 0 ) , TEMP 

I  l=NA-NB+  1 

DO  6  1  Jl  =  1, I  1 

Q  (  Jl)  =0.  ODO 

KKMAX  =  NA-NB+  1 

DO  391  KK= 1 , KKMAX 

K=KK- 1 

TEMP=0.  ODO 

IF(K-l)  301,211,211 

DO  29  1  JJ=1,K 

J=JJ- 1 

I  1=NB-K+J 

IF  (  1 1 )  29  1,22  1,221 
I  2=NA-NB- J 

TEMP=TEMP+B(  I1+1)*Q(I2+1) 

CONTI  NUE 
I  1=NA-NB-K 
I  2  =NA“  K 

Q(  1 1+1)=A ( 12+1) -TEMP 

RETURN 

END 

SUBROUTINE  POL Y ( N , A , R , C , P R , PC , RHO , PH I ) 

DOUBLE  PRECISION  A ( 10 ) , R , C , PR , PC , RHO , PHI , V 1 , V2 , Wl , W2 , T 1 
I F ( RHO )  10,5,10 

R  =  A  (  l) 

C=  0. ODO 
PR=A (2) 

PC=  0. 0  DO 
RETURN 
Vl=  1.0D0 
V2  =  0 . ODO 
R=A  (  l) 

C=  0.  ODO 
PR=0.  ODO 
PC=0.  ODO 

W l=RHO*DCOS (PHI) 

W2=RHO*DSIN (PHI) 
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NN=N+1 

DO  20  1=2, NN 
Tl=Wl*Vl-W2*V2 
V2=W2*Vl+Wl*V2 
VI  =T  1 

R=R+A  ( I ) *V 1 
C=C+A ( I )*V2 
PR=  PR+A ( I )* { I  -  l)*Vl 
20  PC=  PC+A  (I)MI"1)*V2 
PR=  PR/R  HO 
PC=PC/RHO 
5001  RETURN 
END 

FUNCTION  ARCTA  ( X,  Y ) 

SUBPROGRAM  TO  COMPUTE  ARCTANGENT 
VALUES  IN  INTERVAL  -PI  TO  PI 
DOUBLE  PRECISION  PI , HPI , X , Y , ARCTA 
PI= 3. 14159265358979324 
HPI =1.57079632679489662 
I F  (  X )  1,2, 3 

1  ARCT A=  DAT AN ( Y/ X ) +PI*  DS IGN (1.0D0,Y) 
RETURN 

3  A  RCT  A=  DATAN  (  Y/X) 

RET  URN 

2  IF ( Y )  4,5,6 

4  ARCT A= -HPI 
RETURN 

5  ARCTA=  0. 0D0 
RETURN 

6  A  RCT  A=  HPI 
RETURN 
END 

SUBROUTINE  STATS  (AM,  SD,  G) 


STATS  CALCULATES  MEAN  AND  STANDARD  DEVIATION 
OF  A  SET  OF  100  NUMBERS 
DOUBLE  PRECISION  G ( 1 0 0 ) , Si , S2 , RM, AM , SD 
COMMON/F  REQ/FONE ,  FTWO 
COMMON/MAG/RMAG1,  RMAG2 
S  1  =  0.  0D0 
S  2  =  0 . 0  DO 
RM= 100. 0D0 
DO  10  J=  1,  1  00 
10  S 1=S 1+G ( J) 

AM=S 1/RM 
DO  20  J=  1,  10  0 

20  S2=S2+(G(J)-AM)*(G(J)-AM> 

S  D=DSQRT (S2/RM) 

RETURN 

END 
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